Load all important libraries

library(dMod)
library(stringr) # Um bequem mit strings zu arbeiten
library(tidyverse) # Viele Funktionen, u.a. für data.frames und ggplot2 für schöne plots
library(magrittr) # der Pipe-operator %>%: z.B: x = a; y=f(x); z=g(y); wird zu z= a %>% f %>% g
library(conveniencefunctions)
library(libSBML)
load("workspace.rda")

Model setup

Read the model

This is basically copied from a libSBML-example

filename <- "/home/denial/Promotion/Projects/methacetin_fitting/model/met13_pkpd_7.xml" # Attention: tilde expansion doesn't work
d <-  readSBML(filename)

m <- SBMLDocument_getModel(d)
# errors   = SBMLDocument_getNumErrors(d);
# SBMLDocument_printErrors(d);


level   = SBase_getLevel  (d);
version = SBase_getVersion(d);

cat("File: ",filename," (Level ",level,", version ",version,")\n");
cat("  model id: ", ifelse(Model_isSetId(m), Model_getId(m) ,"(empty)"),"\n");

cat( "functionDefinitions: ", Model_getNumFunctionDefinitions(m) ,"\n" );
cat( "    unitDefinitions: ", Model_getNumUnitDefinitions    (m) ,"\n" );
cat( "   compartmentTypes: ", Model_getNumCompartmentTypes   (m) ,"\n" );
cat( "        specieTypes: ", Model_getNumSpeciesTypes       (m) ,"\n" );
cat( "       compartments: ", Model_getNumCompartments       (m) ,"\n" );
cat( "            species: ", Model_getNumSpecies            (m) ,"\n" );
cat( "         parameters: ", Model_getNumParameters         (m) ,"\n" );
cat( " initialAssignments: ", Model_getNumInitialAssignments (m) ,"\n" );
cat( "              rules: ", Model_getNumRules              (m) ,"\n" );
cat( "        constraints: ", Model_getNumConstraints        (m) ,"\n" );
cat( "          reactions: ", Model_getNumReactions          (m) ,"\n" );
cat( "             events: ", Model_getNumEvents             (m) ,"\n" );
cat( "\n" );

Einschub: der Pipe-Operator %>%

# Mit dem Pipe-Operator kann man Funktionen verketten
# Standardmäßig wird das vorherige Ergebnis als erstes Argument von der nächsten Funktion eingesetzt. wenn man das nicht will, kann man es als . woanders hinsetzen

f <- function(x) x^2;
g <- function(x,y) x-y;

2 %>% f # f(2)
2 %>% f %>% g(3) # g(f(2),3) = 4-3 = 1
2 %>% f %>% g(3,.) # g(3,f(2)) = 3-4 = -1

# Man kann auch Funktionen definieren, die mit . losgehen, was dann das Argument für die Funktion ist.
h <- . %>% sqrt %>% add(5)
h(4)
# Das ist besonders nützlich in lapply, sapply und so weiter, wo man über eine liste(oder einen vektor) immer wieder die gleiche funktion laufen lässt (syntactic sugar für for-loops)
sapply(1:4, function(i) i^2 +5)
sapply(1:4 , . %>% raise_to_power(2) %>% add(5) ) # dassselbe

Assignment Rules

  1. Get the assignment rules as string
  2. Apply the assignment rules onto themselves with str_replace. This is for assigments that convert parameters into other parameters
# get rules
nrules <- Model_getNumRules(m)
lrules <- Model_getListOfRules(m)
rules <- structure(sapply(0:(nrules-1), . %>% ListOfRules_get(lrules,.) %>% Rule_getFormula), 
                   names = sapply(0:(nrules-1), . %>% ListOfRules_get(lrules,.) %>% Rule_getId))
rulenames <- names(rules)

# "Cure" rules: Since I do a parameter trafo for the units, I don't want to have any unit conversions via some rules
# A "bad" rule would be eg "QC = CO*3600/100", since I take care of the Units later. Therefore, the rule should be only "QC = CO" 
rules <- rules %>% str_replace_all(c("1000" = "1", "3600" = "1", "\\b60\\b" = 1)) %>% set_names(rulenames)


# Apply the rules onto themselves to insert parameter transformations
# Final goal is to have a named vector where 
  # names are the "inner" parameters that are used within the model
  # values are functions of "outer" parameters that are fed into the model

# apply rules 1st time
rules <- paste0("(", rules, ")") %>% set_names(paste0("\\b", rulenames, "\\b"))
rules <- str_replace_all(rules, rules) %>% set_names(rulenames)

# apply rules 2nd time
rules <- paste0("(", rules, ")") %>% set_names(paste0("\\b", rulenames, "\\b"))
rules <- str_replace_all(rules, rules) %>% set_names(rulenames)

# check if any of the rules are functions of other rules
indices <- rules %>% sapply(. %>% str_detect(paste0("\\b",rulenames, "\\b")) %>% any)
rules %>% extract(indices) %>% sapply(. %>% str_detect(paste0("\\b",rulenames, "\\b")) %>% extract(rulenames,.)) # check works, none of the 

# print(rules)

# getSymbols(rules) # These are the "outer" parameters

Reactions

# get reactions
nreactions <- Model_getNumReactions(m)
lreactions <- Model_getListOfReactions(m)
reactions <- lapply(0:(nreactions-1), function(i) {
  
  reaction <-  ListOfReactions_get(lreactions,i)

  nreactants <- reaction %>% Reaction_getNumReactants()
  if (nreactants > 0) {
    lreactants <- reaction %>% Reaction_getListOfReactants()
    myreactant_stoichiometries <- lapply(0:(nreactants-1), . %>% ListOfSpeciesReferences_get(lreactants,.) %>% SpeciesReference_getStoichiometry)
    myreactant_IDs <- lapply(0:(nreactants-1), . %>% ListOfSpeciesReferences_get(lreactants,.) %>% SimpleSpeciesReference_getSpecies)
    from = paste0(paste0("(", myreactant_stoichiometries, "*", myreactant_IDs, ")"), collapse = "+")
  } else {
    from = ""
  }

  nproducts <- reaction %>% Reaction_getNumProducts()
  if (nproducts > 0) {
    lproducts <- reaction %>% Reaction_getListOfProducts()
    myproduct_stoichiometries <- lapply(0:(nproducts-1), . %>% ListOfSpeciesReferences_get(lproducts,.) %>% SpeciesReference_getStoichiometry)
    myproduct_IDs <- lapply(0:(nproducts-1), . %>% ListOfSpeciesReferences_get(lproducts,.) %>% SimpleSpeciesReference_getSpecies)
    to = paste0(paste0("(", myproduct_stoichiometries, "*", myproduct_IDs, ")"), collapse = "+")
  } else {
    to = ""
  }

  # Apply rules to rate expressions
  myrules <- rules
  mynames <- names(myrules)
  absorption_indices <- str_detect(myrules, "Absorption") 
  absorption_rules <- myrules[absorption_indices] %>% structure(names = mynames[absorption_indices])
  myrules <- structure(myrules[!str_detect(myrules, "Absorption")], names = paste0("\\b", mynames[!str_detect(myrules, "Absorption")], "\\b"))

  rate <-  reaction %>% Reaction_getKineticLaw() %>% KineticLaw_getFormula() %>% str_replace_all(myrules)
  
  description <- reaction %>% Reaction_getId()

  # Incorporate the absorption:
  # For this I add another reaction which absorbs the oral dose, e.g. D_apap_sul
  if(description %>% str_detect("Absorption")) {
    from[2] <- names(absorption_rules)[absorption_rules %>% str_detect(paste0(description, "\\)*$"))]
    to[2] <- ""
    rate[2] <- rate
    description[2] <- description
  }
  # print(i)
  
  return(data.frame(from = from, to  = to, rate = rate, description = description, stringsAsFactors = F))
}) %>% do.call(rbind,.) # make one big data.frame out of it


# Build the dMod-object "eqnlist", which stores the reactions and stoichiometries
el <- eqnlist()
for(i in 1:nrow(reactions)) el <- addReaction(el, reactions$from[i], reactions$to[i], reactions$rate[i], reactions$description[i])


# Convert to "eqnvec", which is basically a named vector of the ODEs and the names denote the states
f <- el %>% as.eqnvec()

# This is the ODE when every rule is applied
print(f %>% head)

Parameters

# get the parameters from the definition
nsbml_parameters <- m %>% Model_getNumParameters()
lsbml_parameters <- m %>% Model_getListOfParameters()
sbml_parameters <- structure(sapply(0:(nsbml_parameters-1), . %>% ListOfParameters_get(lsbml_parameters,.) %>% Parameter_getValue), 
                             names =  sapply(0:(nsbml_parameters-1), . %>% ListOfParameters_get(lsbml_parameters,.) %>% Parameter_getId))
sbml_parameter_descriptions <- structure(sapply(0:(nsbml_parameters-1), . %>% ListOfParameters_get(lsbml_parameters,.) %>% Parameter_getName), 
                                         names = sapply(0:(nsbml_parameters-1), . %>% ListOfParameters_get(lsbml_parameters,.) %>% Parameter_getId)) %>% {.[order(names(.))]}

# unit conversion
sbml_parameter_units <- structure(sapply(0:(nsbml_parameters-1), . %>% ListOfParameters_get(lsbml_parameters,.) %>% Parameter_getUnits), names = sapply(0:(nsbml_parameters-1), . %>% ListOfParameters_get(lsbml_parameters,.) %>% Parameter_getId))
# The factors to bring each parameter to the units seconds, grams, litres
multiplication_factors <- sapply(0:(nsbml_parameters-1), . %>% 
         ListOfParameters_get(lsbml_parameters,.) %>% 
         Parameter_getDerivedUnitDefinition %>%
         {myunitdef <- .
           nunits <- UnitDefinition_getNumUnits(myunitdef)
         lunits <- UnitDefinition_getListOfUnits(myunitdef)
         lapply(0:(nunits-1), . %>%
                  ListOfUnits_get(lunits,.) %>%
                  {(Unit_getMultiplier(.) * 10^(Unit_getScale(.)))^Unit_getExponent(.)}) %>% Reduce("*",.)
         } %>%
         {.}
       )
sbml_parameters <- sbml_parameters * multiplication_factors


# Initial assignments
niass <- m %>% Model_getNumInitialAssignments()
liass <- m %>% Model_getListOfInitialAssignments()
iass <- structure(sapply(0:(niass-1), . %>% ListOfInitialAssignments_get(liass,.) %>% print %>%  str_split("\n") %>% sapply( . %>% {.[str_detect(.,"<ci>")]}) %>% str_replace_all("^.*<ci> ", "") %>% str_replace_all(" </ci>.*$", "")), names = sapply(0:(niass-1), . %>% ListOfInitialAssignments_get(liass,.) %>% InitialAssignment_getId)) 

nspecies <- Model_getNumSpecies(m)
lspecies <- Model_getListOfSpecies(m)
species <- lapply(0:(nspecies-1), . %>% {ListOfSpecies_get(lspecies,.)} %>% Species_getId()) %>% do.call(c,.)

inits <- structure(sapply(0:(nspecies-1), . %>% ListOfSpecies_get(lspecies,.) %>% Species_getInitialAmount), names = species)
inits[names(iass)] <- sbml_parameters[iass]

# All parameters combined, more than actually needed, because many of them have been replaced due to Rules
all_pars <- c(sbml_parameters, inits) # all possible parameters

# To check if the unit conversion worked
data.frame(par = names(sbml_parameters), value = sbml_parameters, original_unit= sbml_parameter_units,  multiplication_factor = multiplication_factors, original_value = sbml_parameters/multiplication_factors, row.names = NULL) %>% 
  head()
#%>%
  # filter(par %in% c("MET2APAP_HLM_CL", "fumic_metc13", "MPPGL", "BW", "FVli"))  #%>% summarise(prod(multiplication_factor))
  # filter(par %in% c("CO", "QC"))
# rules[str_detect(rules, "MET2APAP")]
# rules[str_detect(rules, "CO")]

C-Code

the odemodel() command takes a while because sensitivity equations are calculated for derivatives and then the whole system is compiled into c-code.

# myodemodel <- odemodel(f, modelname = "methacetin") # this compiles the ode into a c-file. you can comment out this and the next line after it has been run once to save time.
# save(myodemodel, file = "methacetin.rda")
load("methacetin.rda")

Prediction function

Make a prediction function from the odemodel. x will be a function x(times, pars)

x <- Xs(myodemodel) # make prediction function
loadDLL(x)

# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]

Example plot

This plot is supposed to be the first plot from chunk 3 of the html-file that you sent me via email.

pars["Ave_metc13"] <- 0
pars["D_metc13"] <- 0.1 #100mg

mytimes = seq(0,8*3600, length.out = 80) #8 hours

# pred <- (g*x)(0:50, mypars, deriv = F, conditions = rownames(attr(data, "cond")))
prediction <- (x)(times = mytimes, pars = pars, deriv = F, conditions = "c1")

# Compute the volumes in litres for each state
volumes <- lapply(0:(nspecies-1), . %>% {ListOfSpecies_get(lspecies,.)} %>% Species_getCompartment()) %>% do.call(c,.) %>% set_names(species)
myrules <- rules %>% set_names(paste0("\\b", rulenames, "\\b"))
volumes <- str_replace_all(volumes, myrules) %>% set_names(species)
vol_fun <- funC0(volumes)
vol <- do.call(vol_fun, as.list(pars)) %>% t %>% {data.frame(volume=., name = rownames(.), stringsAsFactors = F)}

# Plot
plot_fun <- function(pred) pred %>% 
  wide2long() %>% 
  full_join(vol,by= "name") %>% 
  mutate(value = value * 1000, time = time/3600, concentration = value/volume) %>% # scale g to mg, s to h
  separate(col = name, into = c("compartment", "substance"), sep = "_", extra = "merge") %>% 
  filter(compartment %in% c("Ali", "Ave", "Alu")) %>% # plot only liver, venuous and lung (as in the html file)
  ggplot(aes(x= time, y = concentration)) +
  geom_line(aes(color = compartment)) + 
  facet_wrap(~substance, scales = "free") + 
  theme_dMod() + scale_color_dMod()



plot_fun(prediction)

prediction_long <- (x)(times = mytimes*20,pars = pars, deriv = F, conditions = "c1")
plot_fun(prediction_long)

# Second plot from chunk 3
pars <- all_pars[getParameters(x)]
pars["Ave_metc13"] <- 0.1
pars["D_metc13"] <- 0 #100mg
# pars["Kplu_metc13"] <- 1
# pars["FVlu"] <- 7.6*10^(-6) * 10

prediction_second_plot <- (x)(times = mytimes,pars = pars, deriv = F, conditions = "c1")
plot_fun(prediction_second_plot)

Modelling

Reduce the model complexity by inserting the fixed parameter values

free_parameters <- c("APAPGLU_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPSUL_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPCYS_HLM_CL",  # Vmax value
                     "APAPCYS_Km"  # Km value
                     )


fixed_parameters <- pars[!(names(pars)%in%c(free_parameters,names(f)[1]))] %>% names


myodemodel <- odemodel(f, modelname = "methacetin_reduced", fixed = fixed_parameters)

x <- Xs(myodemodel)

Observation function

observables <- c(apap = "Ave_apap/(BW*FVve)", 
                 apap_glu = "Ave_apap_glu/(BW*FVve)", 
                 apap_sul = "Ave_apap_sul/(BW*FVve)", 
                 apap_cys = "Ave_apap_cys/(BW*FVve)"
                 )

g <- Y(observables, x, parameters = free_parameters)
(g*x)(mytimes, pars) %>% set_names("condition1") %>% plotPrediction(name %in% names(observables))

Data

Load the data and transform it

myfiles <- list.files("~/Promotion/Projects/methacetin_fitting/data/", full.names = T)


raw_data <- myfiles %>% lapply(. %>% read.table(header = T, sep = "\t", stringsAsFactors = F))


data <-
  raw_data %>% 
  lapply(. %>% 
           select(-contains("_mol")) %>% 
           gather("name_std", "std", ends_with("_sd")) %>% 
           mutate(name_std = str_replace(name_std, "_sd","")) %>% 
           gather("name_sigma", "sigma", ends_with("_se")) %>% 
           mutate(name_sigma = str_replace(name_sigma, "_se","")) %>% 
           {gather(.,"name", "value", one_of(.$name_std))} %>% 
           filter(name == name_std, name == name_sigma) %>% 
           # select(name,time,value,sigma) %>% 
           {.}) %>% 
    do.call(dMod::combine,.) %>% 
    mutate(D_apap = "D_apap", Ave_apap = "Ave_apap" ) %>% 
    {.$D_apap[.$study=="Chan1997"] <- 1400 / 1000
    .$Ave_apap[.$study=="Chan1997"] <- 0
    .$D_apap[.$study=="Chiew2010"] <- 5600 / 1000
    .$Ave_apap[.$study=="Chiew2010"] <- 0
    .$D_apap[.$study=="Critchley2005"] <- 1400 /1000
    .$Ave_apap[.$study=="Critchley2005"] <- 0
    .$D_apap[.$study=="Rawlins1977"] <- .$dose[.$study=="Rawlins1977"] * (.$route[.$study=="Rawlins1977"]%in%"oral") / 1000
    .$Ave_apap[.$study=="Rawlins1977"] <- .$dose[.$study=="Rawlins1977"] * (.$route[.$study=="Rawlins1977"]%in%"iv") / 1000
    .
    } %>% 
  mutate(time = time * 3600, value = value/1000, sigma = sigma/1000) %>% 
  select(-group, -health_status, - name_std, - name_sigma, -std, -ethnicity, -route, -dose, -substance) %>%
  # filter(!is.na(sigma)) %>% 
  # as.datalist() %>% 
           {.}
mydatalist <- data %>% filter(!is.na(sigma)) %>% select(-n) %>% as.datalist()

plot(mydatalist)

Parameter transformations to define the conditions

conditions <- mydatalist %>% attr("condition.grid")

p_list <- lapply(1:nrow(conditions), function(i) {
  trafo <- as.character(pars) %>% set_names(names(pars))
  cond <- unlist(conditions[i,])[2:3]
  
  trafo[names(cond)] <- cond
  trafo[free_parameters] <- paste0("exp(log", free_parameters, ")")
  
  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p <- NULL
for(i in 1:length(p_list)) { p <<- p + p_list[[i]]}

pouter <- log(pars[free_parameters]) %>% set_names(paste0("log",names(.)))


mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
plotCombined(mypred, mydatalist, name%in% names(observables))

Fitting

obj <- normL2(mydatalist, (g*x*p))


# myfit <- mstrust(objfun = obj, center = pouter, studyname = "methacetin", cores = 3)
# save(myfit, file = "fit.rda")
load("fit.rda")
fitted_pars <- myfit %>% as.parframe() %>% as.parvec()
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
mypred <- (g*x*p)(seq(0, 48*3600, length.out = 100), fitted_pars)

liver <- names(f)[str_detect(names(f), "li")]
medication <- c(names(f)[str_detect(names(f), "D_apap$")], "Ave_apap")

myfit %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
myfit %>% as.parframe() %>% plotValues()

plotCombined(mypred, mydatalist, name %in% c(names(observables), liver, medication))

##   index    value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1     8 2407.821      TRUE         60         -5.125968     -2.129845
## 2    35 2407.821      TRUE         71         -5.125973     -2.129851
## 3    21 2407.821      TRUE         39         -5.125965     -2.129842
## 4    29 2407.821      TRUE         67         -5.125962     -2.129839
## 5    33 2407.821      TRUE         76         -5.125966     -2.129843
## 6    22 2407.821      TRUE         75         -5.125969     -2.129847
##   logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km
## 1         -6.118641         -10.49793     -6.301154
## 2         -6.118641         -10.49793     -6.301161
## 3         -6.118641         -10.49794     -6.301223
## 4         -6.118641         -10.49794     -6.301221
## 5         -6.118641         -10.49794     -6.301214
## 6         -6.118641         -10.49794     -6.301198
## logAPAPCYS_HLM_CL:  2.75936204923354e-05
##     logAPAPCYS_Km:  0.00183418625326802
## logAPAPGLU_HLM_CL:  0.00594046463001759
##     logAPAPGLU_Km:  0.118855661904664
## logAPAPSUL_HLM_CL:  0.00220144570190944

Free some additional parameters

Free other parameters try 1 - not good

Here I freed some additional parameters, but since I freed “Ka_apap_sul” instead of “Ka_apap”, the results don’t make much additional sense.

load("methacetin.rda")

x <- Xs(myodemodel) # make prediction function
loadDLL(x)

# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]

free_parameters1 <- c("APAPGLU_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPSUL_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPCYS_HLM_CL",  # Vmax value
                     "APAPCYS_Km",  # Km value

                     "Ka_apap_sul", "F_apap_sul",
                     "Kpre_apap", "Kpki_apap", "Kpli_apap",
                     "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
                     "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
                     "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu",
                     "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
                     "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
                     )

fixed_parameters1 <- pars[!(names(pars)%in%c(free_parameters1,names(f)[1]))] %>% names

mydatalist <- data  %>% filter(!is.na(sigma)) %>% as.datalist()

conditions <- mydatalist %>% attr("condition.grid")

p_list <- lapply(1:nrow(conditions), function(i) {
  trafo <- as.character(pars) %>% set_names(names(pars))
  cond <- unlist(conditions[i,])[2:3]

  trafo[names(cond)] <- cond
  trafo[free_parameters1] <- paste0("exp(log", free_parameters1, ")")

  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p1 <- NULL
for(i in 1:length(p_list)) { p1 <<- p1 + p_list[[i]]}

pouter <- log(pars[free_parameters1]) %>% set_names(paste0("log",names(.)))


# mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
# plotCombined(mypred, mydatalist, name%in% names(observables))

obj1 <- normL2(mydatalist, (g*x*p1))

# job1 <- runbg({myfit <- mstrust(objfun = obj1, center = pouter, studyname = "methacetin", cores = 20, fits = 100); myfit}, machine = "knecht3", filename = "job1")
# myfit1 <- job1$get()$knecht3
# save(myfit1, file = "myfit1.rda")
# job1$purge()
load("myfit1.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit1 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit1 %>% as.parframe())+scale_y_log10()

mypred1 <- (g*x*p1)(mytimes*4, myfit1 %>% as.parframe() %>% as.parvec)
plotCombined(mypred1, mydatalist, name %in% names(observables))

##   index    value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1    71 2749.704     FALSE        100         -3.121425     0.9276677
## 2    17 2871.109     FALSE        100         -4.825791    -0.8966849
## 3     4 2879.499      TRUE         82         -4.836514    -2.0500243
## 4    83 3006.360      TRUE         38         -4.810375    -1.6071039
## 5    11 3025.134      TRUE         47         -5.068367    -2.1571283
## 6    60 3063.915     FALSE        100         -4.778175    -1.9189665
##   logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap_sul
## 1         -5.966762         -8.877054     -8.056205      -9.043602
## 2         -7.312292         -9.776533     -6.317704      -7.411144
## 3         -6.424722         -9.745092     -7.247758      -7.449345
## 4         -6.432805        -10.767499     -6.451981      -6.734894
## 5         -6.368207         -9.236633     -7.298571      -8.975554
## 6         -6.479809        -10.151009     -5.692777      -5.931667
##   logF_apap_sul logKpre_apap logKpki_apap logKpli_apap logKpre_apap_cys
## 1    0.53130523   -0.2895295   -0.9819923   0.32392147        3.4721895
## 2    0.07956807   -0.3510074    0.2160174   0.99467641        0.7860616
## 3   -0.52122381   -0.3932631   -1.0615523  -0.09235355        0.9514369
## 4    0.16252315   -0.3785620    0.1383524   0.27175050       -0.8202826
## 5    1.21332440   -0.2100422   -1.7021606  -0.01216570        1.9565013
## 6   -2.38273975   -0.4474130   -1.6230284   0.04538435        0.4837584
##   logKpki_apap_cys logKpli_apap_cys logKpre_apap_glu logKpki_apap_glu
## 1       -0.9981306       -0.7078988       -1.1713447       -0.4586217
## 2        1.4876644        0.5391938        0.5938409        0.6894229
## 3        1.3055012        0.4413984        0.7986194        0.3810163
## 4       -0.1222928       -0.2784549       -0.6756876        0.1384249
## 5        1.2069110        0.6198618        0.5661661       -1.1336736
## 6       -0.1000206       -0.1492559        0.6852967       -0.2904269
##   logKpli_apap_glu logKpre_apap_sul logKpre_co2c13 logKpli_co2c13
## 1        1.2820601        0.8610661     -1.8050599     -1.0958784
## 2       -1.2743622       -2.0279402     -1.6334736      0.2964431
## 3       -0.7746656       -1.5237195      0.3130873      0.3634344
## 4        0.4949366       -0.5271704     -0.8464643     -1.1150170
## 5       -1.4231745       -1.4809556      0.9901271     -0.4769397
## 6        1.7057846       -1.3223539     -1.2621168     -0.9185781
##   logKpre_metc13 logKpli_metc13
## 1     -2.4242875     -0.4355059
## 2     -1.3074538     -0.5878608
## 3      0.0672823      1.2834798
## 4     -2.5253310      1.1213766
## 5     -1.0586655      1.3624494
## 6     -1.3302118      1.4817382
## logAPAPCYS_HLM_CL:  0.000139554700132717
##     logAPAPCYS_Km:  0.000317128135451512
## logAPAPGLU_HLM_CL:  0.0440942903336689
##     logAPAPGLU_Km:  2.52860494207782
## logAPAPSUL_HLM_CL:  0.00256252477737589
##     logF_apap_sul:  1.70115125864559
##    logKa_apap_sul:  0.000118144550067709
##      logKpki_apap:  0.374564129323335
##  logKpki_apap_cys:  0.368567801198445
##  logKpki_apap_glu:  0.63215432668864
##      logKpli_apap:  1.38253873675798
##  logKpli_apap_cys:  0.492678343277563
##  logKpli_apap_glu:  3.60405672997012
##    logKpli_co2c13:  0.334245860035197
##    logKpli_metc13:  0.646937302760331
##      logKpre_apap:  0.748615739138366
##  logKpre_apap_cys:  32.2071824118807
##  logKpre_apap_glu:  0.309949884477503
##  logKpre_apap_sul:  2.36568138653131
##    logKpre_co2c13:  0.164464609213266
##    logKpre_metc13:  0.0885411804880305

Free other parameters 2 - not good

Here I freed some less parameters than in try 1 but it suffers from the same problem

load("methacetin.rda")

x <- Xs(myodemodel) # make prediction function
loadDLL(x)

# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]

free_parameters2 <- c("APAPGLU_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPSUL_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPCYS_HLM_CL",  # Vmax value
                     "APAPCYS_Km",  # Km value

                     "Ka_apap_sul", "F_apap_sul"    #,   # total dumm, Ka und F sind so redundant wie es nur geht.
                     # "Kpre_apap", "Kpki_apap", "Kpli_apap",
                     # "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
                     # "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
                     # "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu",
                     # "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
                     # "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
                     )

fixed_parameters2 <- pars[!(names(pars)%in%c(free_parameters2,names(f)[1]))] %>% names

mydatalist <- data  %>% filter(!is.na(sigma)) %>% as.datalist()

conditions <- mydatalist %>% attr("condition.grid")

p_list <- lapply(1:nrow(conditions), function(i) {
  trafo <- as.character(pars) %>% set_names(names(pars))
  cond <- unlist(conditions[i,])[2:3]

  trafo[names(cond)] <- cond
  trafo[free_parameters2] <- paste0("exp(log", free_parameters2, ")")

  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p2 <- NULL
for(i in 1:length(p_list)) { p2 <<- p2 + p_list[[i]]}

pouter <- log(pars[free_parameters2]) %>% set_names(paste0("log",names(.)))


# mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
# plotCombined(mypred, mydatalist, name%in% names(observables))

obj2 <- normL2(mydatalist, (g*x*p2))

# job2 <- runbg({myfit <- mstrust(objfun = obj2, center = pouter, studyname = "methacetin", cores = 12, fits = 100); myfit}, machine = "knecht4", filename = "job2")
# 
# myfit2 <- job2$get()$knecht4
# save(myfit2, file = "myfit2.rda")
# job2$purge()
load("myfit2.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit2 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit2 %>% as.parframe())+scale_y_log10()

mypred2 <- (g*x*p2)(mytimes*4, myfit2 %>% as.parframe() %>% as.parvec)
plotCombined(mypred2, mydatalist, name %in% names(observables))

##   index    value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1    24 3134.039      TRUE         25         -5.067594     -2.064909
## 2    56 3134.040      TRUE         30         -5.067589     -2.064904
## 3    20 3134.040      TRUE         26         -5.067595     -2.064908
## 4    89 3134.040      TRUE         16         -5.067592     -2.064905
## 5    37 3134.040      TRUE         86         -5.067589     -2.064901
## 6    32 3134.040      TRUE         20         -5.067590     -2.064902
##   logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap_sul
## 1         -6.117757         -10.48519     -6.241245      -7.346186
## 2         -6.117757         -10.48520     -6.241296      -7.021465
## 3         -6.117758         -10.48501     -6.240447      -7.216762
## 4         -6.117758         -10.48501     -6.240447      -6.306130
## 5         -6.117757         -10.48501     -6.240481      -7.986534
## 6         -6.117757         -10.48501     -6.240461      -7.899805
##   logF_apap_sul
## 1     41.111322
## 2     43.933670
## 3     50.970546
## 4     43.726021
## 5      7.397659
## 6     34.039840
## logAPAPCYS_HLM_CL:  2.79472212506865e-05
##     logAPAPCYS_Km:  0.00194743024961977
## logAPAPGLU_HLM_CL:  0.00629755070545838
##     logAPAPGLU_Km:  0.126829778557011
## logAPAPSUL_HLM_CL:  0.00220339319347716
##     logF_apap_sul:  715187883400265600
##    logKa_apap_sul:  0.000645048206917656

Free other parameters 3 - not good

Same problem as in tries 1 and 2

load("methacetin.rda")

x <- Xs(myodemodel) # make prediction function
loadDLL(x)

# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]

free_parameters3 <- c("APAPGLU_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPSUL_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPCYS_HLM_CL",  # Vmax value
                     "APAPCYS_Km",  # Km value

                     "Ka_apap_sul", "F_apap_sul"    ,
                     "Kpre_apap", "Kpki_apap", "Kpli_apap",
                     "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
                     "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
                     "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
                     # "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
                     # "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
                     )

fixed_parameters3 <- pars[!(names(pars)%in%c(free_parameters3,names(f)[1]))] %>% names

mydatalist <- data  %>% filter(!is.na(sigma)) %>% as.datalist()

conditions <- mydatalist %>% attr("condition.grid")

p_list <- lapply(1:nrow(conditions), function(i) {
  trafo <- as.character(pars) %>% set_names(names(pars))
  cond <- unlist(conditions[i,])[2:3]

  trafo[names(cond)] <- cond
  trafo[free_parameters3] <- paste0("exp(log", free_parameters3, ")")

  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p3 <- NULL
for(i in 1:length(p_list)) { p3 <<- p3 + p_list[[i]]}

pouter <- log(pars[free_parameters3]) %>% set_names(paste0("log",names(.)))


# mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
# plotCombined(mypred, mydatalist, name%in% names(observables))

obj3 <- normL2(mydatalist, (g*x*p3))

# job3 <- runbg({myfit <- mstrust(objfun = obj3, center = pouter, studyname = "methacetin", cores = 15, fits = 100); myfit}, machine = "knecht1", filename = "job3")
# myfit3 <- job3$get()$knecht1
# save(myfit3, file = "myfit3.rda")
# job3$purge()
load("myfit3.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit3 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit3 %>% as.parframe())+scale_y_log10()

mypred3 <- (g*x*p3)(mytimes*4, myfit3 %>% as.parframe() %>% as.parvec)
plotCombined(mypred3, mydatalist, name %in% names(observables))

##   index    value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1    59 2562.482      TRUE         84         -5.201098     -2.097660
## 2    60 2689.397      TRUE         48         -5.383774     -2.043510
## 3    55 2856.425     FALSE        100         -5.109883     -1.468816
## 4    48 3019.107     FALSE        100         -6.865873     -5.295360
## 5     1 3039.058     FALSE        100         -4.638935     -1.383292
## 6    89 3066.281      TRUE         76         -6.464966     -4.620926
##   logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap_sul
## 1         -6.555927        -10.709675   -6.33769669      -6.042436
## 2         -6.862922        -10.910624   -6.20766738      -8.289160
## 3         -6.656348         -9.872593   -4.75970985      -6.579562
## 4         -6.050377         -9.618928   -6.86802281      -6.542891
## 5         -6.709802         -7.521883   -0.09424782      -7.980098
## 6         -5.857182        -10.682944   -8.86895808      -7.943206
##   logF_apap_sul logKpre_apap logKpki_apap logKpli_apap logKpre_apap_cys
## 1   2.808457812   -0.3669234    0.7613894    0.1672635       -0.6048195
## 2  -1.512335536   -0.1336978   -0.2052586    0.5449920       -1.5425344
## 3   0.009716211   -0.2655233   -1.6898116    0.6768234       -0.2769801
## 4  -0.661252544   -0.1334166    0.2229806   -0.6106942       -0.4045914
## 5   2.748915536   -0.4468966   -0.6003575    0.3837791       -0.3193860
## 6   1.679453017   -0.1192626   -0.5590351   -0.6042076       -2.9191368
##   logKpki_apap_cys logKpli_apap_cys logKpre_apap_glu logKpki_apap_glu
## 1       -0.1241948      -0.04614371       -0.8681152       0.36911867
## 2       -0.2511150       1.52077842       -0.8083817       0.16712529
## 3        0.7221217      -0.09379088       -1.0941684       0.26637337
## 4        2.2752818       0.30050894       -1.5193795       1.64286021
## 5       -2.5272957       1.04456393        0.3086699       0.04221941
## 6        0.4521719      -0.72613317        0.1994445       0.85658306
##   logKpli_apap_glu logKpre_apap_sul
## 1     -0.003140001       -1.0894804
## 2     -0.342208285       -2.2243255
## 3     -0.207101279       -0.1223363
## 4     -0.702140533       -2.0416762
## 5      0.408290615       -1.1143709
## 6     -0.820070493       -0.7383285
## logAPAPCYS_HLM_CL:  2.23278620657747e-05
##     logAPAPCYS_Km:  0.00176837065068179
## logAPAPGLU_HLM_CL:  0.00551051224903083
##     logAPAPGLU_Km:  0.122743302983146
## logAPAPSUL_HLM_CL:  0.0014216640114061
##     logF_apap_sul:  16.5843223443174
##    logKa_apap_sul:  0.00237576415394654
##      logKpki_apap:  2.14124916301601
##  logKpki_apap_cys:  0.883207793863882
##  logKpki_apap_glu:  1.44645924035911
##      logKpli_apap:  1.18206573439193
##  logKpli_apap_cys:  0.954904720644904
##  logKpli_apap_glu:  0.996864923604179
##      logKpre_apap:  0.692862701926673
##  logKpre_apap_cys:  0.546172985593896
##  logKpre_apap_glu:  0.419741938691713
##  logKpre_apap_sul:  0.336391236317653

Free other parameters 4 - not good

Same problem as in tries 1,2,3. The only difference to try 2 is that I removed the structural non-identifiability of Ka_apap_sul and F_apap_sul.

load("methacetin.rda")

x <- Xs(myodemodel) # make prediction function
loadDLL(x)

# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]

free_parameters4 <- c("APAPGLU_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPSUL_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPCYS_HLM_CL",  # Vmax value
                     "APAPCYS_Km",  # Km value

                     "Ka_apap_sul"#, "F_apap_sul"    ,
                     # "Kpre_apap", "Kpki_apap", "Kpli_apap",
                     # "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
                     # "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
                     # "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
                     # "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
                     # "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
                     )

fixed_parameters4 <- pars[!(names(pars)%in%c(free_parameters4,names(f)[1]))] %>% names

mydatalist <- data  %>% filter(!is.na(sigma)) %>% as.datalist()

conditions <- mydatalist %>% attr("condition.grid")

p_list <- lapply(1:nrow(conditions), function(i) {
  trafo <- as.character(pars) %>% set_names(names(pars))
  cond <- unlist(conditions[i,])[2:3]

  trafo[names(cond)] <- cond
  trafo[free_parameters4] <- paste0("exp(log", free_parameters4, ")")

  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p4 <- NULL
for(i in 1:length(p_list)) { p4 <<- p4 + p_list[[i]]}

pouter <- log(pars[free_parameters4]) %>% set_names(paste0("log",names(.)))


# mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
# plotCombined(mypred, mydatalist, name%in% names(observables))

obj4 <- normL2(mydatalist, (g*x*p4))

# job4 <- runbg({myfit <- mstrust(objfun = obj4, center = pouter, studyname = "methacetin", cores = 15, fits = 100); myfit}, machine = "knecht1", filename = "job4")

# save(job4, file = "job4.rda")
# myfit4 <- job4$get()$knecht1
# save(myfit4, file = "myfit4.rda")
# job4$purge()
load("myfit4.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit4 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit4 %>% as.parframe())+scale_y_log10()

mypred4 <- (g*x*p4)(mytimes*4, myfit4 %>% as.parframe() %>% as.parvec, deriv = F)
myplot4 <- plotCombined(mypred4, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot4)
# myplot4
##   index    value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1    65 3134.039      TRUE         31         -5.067596     -2.064909
## 2    47 3134.040      TRUE         74         -5.067589     -2.064901
## 3    51 3134.040      TRUE         82         -5.067591     -2.064903
## 4    49 3134.040      TRUE         58         -5.067589     -2.064902
## 5    17 3134.040      TRUE         18         -5.067592     -2.064904
## 6    91 3134.040      TRUE         25         -5.067590     -2.064903
##   logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap_sul
## 1         -6.117758         -10.48501     -6.240449     69.0940136
## 2         -6.117757         -10.48502     -6.240488     -1.9728623
## 3         -6.117757         -10.48501     -6.240450     -1.9723416
## 4         -6.117757         -10.48501     -6.240446     -0.7434353
## 5         -6.117757         -10.48501     -6.240455     21.5604671
## 6         -6.117758         -10.48499     -6.240395     39.2232311
## logAPAPCYS_HLM_CL:  2.79524088116887e-05
##     logAPAPCYS_Km:  0.00194898079879723
## logAPAPGLU_HLM_CL:  0.0062975431281648
##     logAPAPGLU_Km:  0.126829823077826
## logAPAPSUL_HLM_CL:  0.00220339066957284
##    logKa_apap_sul:  1.01659705294599e+30

Further hypotheses

Freeing the parameters in the upper four fits didn’t help much. On the other hand, lots of the fits didn’t converge. From here, I could do several things.

  1. Introduce scaling factors for each condition
  2. Fill the empty sigmas with fitted sigmas
  3. Run a profile likelihood analysis for the fits 1-4
  4. Re-run the fits 2 and 3 without the structural non-identifiability of Ka_apap_sul and F_apap_sul
  5. If I re-run the fits, I can take the fitted values of the original free parameters as center for the sampling
  6. In data, drop the column n before converting it to a datalist
  7. In the upper fits, I freed the parameters Ka_apap_sul, which of course doesn’t make sense, since no apap_sul is given. I have to free the parameters Ka_apap

Free other parameters 5 - not good

Here I freed Ka_apap which makes sense conceptually (increase absorption rate). But it doesn’t improve the fits much.

load("methacetin.rda")

x <- Xs(myodemodel) # make prediction function
loadDLL(x)

# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]

free_parameters5 <- c("APAPGLU_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPSUL_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPCYS_HLM_CL",  # Vmax value
                     "APAPCYS_Km",  # Km value

                     "Ka_apap"#, "F_apap_sul"    ,
                     # "Kpre_apap", "Kpki_apap", "Kpli_apap",
                     # "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
                     # "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
                     # "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
                     # "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
                     # "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
                     )

fixed_parameters5 <- pars[!(names(pars)%in%c(free_parameters5,names(f)[1]))] %>% names

mydatalist <- data  %>% filter(!is.na(sigma)) %>% as.datalist()

conditions <- mydatalist %>% attr("condition.grid")

p_list <- lapply(1:nrow(conditions), function(i) {
  trafo <- as.character(pars) %>% set_names(names(pars))
  cond <- unlist(conditions[i,])[2:3]

  trafo[names(cond)] <- cond
  trafo[free_parameters5] <- paste0("exp(log", free_parameters5, ")")

  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p5 <- NULL
for(i in 1:length(p_list)) { p5 <<- p5 + p_list[[i]]}

pouter <- log(pars[free_parameters5]) %>% set_names(paste0("log",names(.)))


# mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
# plotCombined(mypred, mydatalist, name%in% names(observables))

obj5 <- normL2(mydatalist, (g*x*p5))

# job5 <- runbg({myfit <- mstrust(objfun = obj5, center = pouter, studyname = "methacetin", cores = 12, fits = 100); myfit}, machine = "knecht1", filename = "job5")

# save(job5, file = "job5.rda")
# myfit5 <- job5$get()$knecht1
# save(myfit5, file = "myfit5.rda")
# job5$purge()
load("myfit5.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit5 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit5 %>% as.parframe())+scale_y_log10()

mypred5 <- (g*x*p5)(mytimes*5, myfit5 %>% as.parframe() %>% as.parvec, deriv = F)
myplot5 <- plotCombined(mypred5, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot5)
# myplot5
##   index    value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1    33 3092.015     FALSE        100         -5.048048     -2.013030
## 2    12 3092.025      TRUE         51         -5.047014     -2.013058
## 3    99 3092.031     FALSE        100         -5.048425     -2.014686
## 4    60 3092.117      TRUE         96         -5.047379     -2.011610
## 5    88 3092.146      TRUE         49         -5.047300     -2.011395
## 6    26 3092.203      TRUE         40         -5.048943     -2.016184
##   logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1         -6.142913         -10.48392     -6.226134  -6.853661
## 2         -6.141722         -10.48302     -6.222827  -6.879080
## 3         -6.141694         -10.48360     -6.225265  -6.880146
## 4         -6.143583         -10.48343     -6.224064  -6.838579
## 5         -6.143710         -10.48332     -6.223574  -6.835689
## 6         -6.140809         -10.48370     -6.225923  -6.898840
## logAPAPCYS_HLM_CL:  2.79827234659552e-05
##     logAPAPCYS_Km:  0.00197708080993727
## logAPAPGLU_HLM_CL:  0.00642185876478594
##     logAPAPGLU_Km:  0.133583353550372
## logAPAPSUL_HLM_CL:  0.0021486564154929
##        logKa_apap:  0.00105558386269278

Free other parameters 6 - not enough fits converged

load("methacetin.rda")

x <- Xs(myodemodel) # make prediction function
loadDLL(x)

# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]

free_parameters6 <- c("APAPGLU_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPSUL_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPCYS_HLM_CL",  # Vmax value
                     "APAPCYS_Km",  # Km value

                     "Ka_apap", #"F_apap_sul"    ,
                     "Kpre_apap", "Kpki_apap", "Kpli_apap",
                     "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
                     "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
                     "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
                     # "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
                     # "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
                     )

fixed_parameters6 <- pars[!(names(pars)%in%c(free_parameters6,names(f)[1]))] %>% names

mydatalist <- data  %>% filter(!is.na(sigma)) %>% as.datalist()

conditions <- mydatalist %>% attr("condition.grid")

p_list <- lapply(1:nrow(conditions), function(i) {
  trafo <- as.character(pars) %>% set_names(names(pars))
  cond <- unlist(conditions[i,])[2:3]

  trafo[names(cond)] <- cond
  trafo[free_parameters6] <- paste0("exp(log", free_parameters6, ")")

  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p6 <- NULL
for(i in 1:length(p_list)) { p6 <<- p6 + p_list[[i]]}

pouter <- log(pars[free_parameters6]) %>% set_names(paste0("log",names(.)))
best_fit <- myfit %>% as.parframe() %>% as.parvec()
pouter[names(best_fit)] <- best_fit

# mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
# plotCombined(mypred, mydatalist, name%in% names(observables))

obj6 <- normL2(mydatalist, (g*x*p6))

# job6 <- runbg({myfit <- mstrust(objfun = obj6, center = pouter, studyname = "methacetin", cores = 10, fits = 100); myfit}, machine = "knecht3", filename = "job6")

# save(job6, file = "job6.rda")
# job6$check()
# myfit6 <- job6$get()$knecht3
# save(myfit6, file = "myfit6.rda")
# job6$purge()
load("myfit6.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit6 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit6 %>% as.parframe())+scale_y_log10()

mypred6 <- (g*x*p6)(mytimes*4, myfit6 %>% as.parframe() %>% {.[2,]} %>% as.parvec)
myplot6 <- plotCombined(mypred6, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot6)
# myplot6
##   index    value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1    16 2654.322      TRUE         88         -6.288191    -4.0064053
## 2    60 3025.841     FALSE        100          5.036281     8.1988600
## 3    39 3059.748     FALSE        100         -3.427671     0.2086894
## 4    19 3210.587      TRUE        100         -4.251426    -1.6417543
## 5    47 3233.777      TRUE         48         -5.809100    -2.7455612
## 6    48 3302.297     FALSE        100         -4.770855    -0.8605885
##   logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1         -6.242811         -9.413691     -6.941120  -5.245524
## 2         -5.281948        -10.686351     -6.408999  -7.197151
## 3         -6.146221         -9.736307     -7.455677  -6.582344
## 4         -4.910387        -10.540902     -7.410143  -7.443229
## 5         -6.315936        -10.030452     -1.701016  -7.739147
## 6         -7.215143         -9.338544     -1.877134  -7.204190
##   logKpre_apap logKpki_apap logKpli_apap logKpre_apap_cys logKpki_apap_cys
## 1  -0.04690208    0.1358953   -0.3757922       2.19032710       -0.4659840
## 2  -0.27191200    1.4639032   -0.5538510      -0.52221741       -0.3912030
## 3  -0.42939025    0.3304541    0.2229900       0.05838160        1.8274448
## 4  -0.18759530   -1.7257917   -0.7659503      -0.40658360        0.1236754
## 5  -0.08418012    0.4028411    0.2491479      -3.17511619       -2.8534696
## 6  -0.34371554   -1.2734140    1.0683678       0.03096511       -1.5696406
##   logKpli_apap_cys logKpre_apap_glu logKpki_apap_glu logKpli_apap_glu
## 1     -0.298471419       -1.2802038       0.54299839       0.09378855
## 2      0.276192230       -1.6506727      -0.45872701      -0.11353331
## 3      0.279366410       -1.3550508      -0.31509066      -0.10951404
## 4      0.009838299       -0.8169902      -0.19239433      -0.32722341
## 5     -0.007563093       -2.0011674       0.04291508      -0.71993332
## 6     -0.671261233        0.5743076      -1.17050260       1.01329858
##   logKpre_apap_sul
## 1      -1.50864561
## 2       0.48599584
## 3      -0.05666556
## 4       0.78554977
## 5      -0.54997601
## 6      -1.66492170
## logAPAPCYS_HLM_CL:  8.15992019166072e-05
##     logAPAPCYS_Km:  0.000967185735063076
## logAPAPGLU_HLM_CL:  0.00185811908010114
##     logAPAPGLU_Km:  0.0181986967501468
## logAPAPSUL_HLM_CL:  0.00194438305778717
##        logKa_apap:  0.0052710585929262
##      logKpki_apap:  1.14556190846099
##  logKpki_apap_cys:  0.627517333300021
##  logKpki_apap_glu:  1.72115983885983
##      logKpli_apap:  0.686745011371109
##  logKpli_apap_cys:  0.741951487277269
##  logKpli_apap_glu:  1.09832748267348
##      logKpre_apap:  0.954180830717175
##  logKpre_apap_cys:  8.93813625877344
##  logKpre_apap_glu:  0.277980630583032
##  logKpre_apap_sul:  0.22120937811377

Free other parameters 6_1 - more iterations didn’t help, it’s still bad

In comparison to try 6, I also drop column n before converting the data to a datalist. I use the best fit from try 6 as center for the sampling

load("methacetin.rda")

x <- Xs(myodemodel) # make prediction function
loadDLL(x)

# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]

free_parameters6 <- c("APAPGLU_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPSUL_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPCYS_HLM_CL",  # Vmax value
                     "APAPCYS_Km",  # Km value

                     "Ka_apap", #"F_apap_sul"    ,
                     "Kpre_apap", "Kpki_apap", "Kpli_apap",
                     "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
                     "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
                     "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
                     # "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
                     # "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
                     )

fixed_parameters6 <- pars[!(names(pars)%in%c(free_parameters6,names(f)[1]))] %>% names

mydatalist <- data  %>% filter(!is.na(sigma)) %>% select(-n) %>% as.datalist()

conditions <- mydatalist %>% attr("condition.grid")

p_list <- lapply(1:nrow(conditions), function(i) {
  trafo <- as.character(pars) %>% set_names(names(pars))
  cond <- unlist(conditions[i,])[2:3]

  trafo[names(cond)] <- cond
  trafo[free_parameters6] <- paste0("exp(log", free_parameters6, ")")

  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p6_1 <- NULL
for(i in 1:length(p_list)) { p6_1 <<- p6_1 + p_list[[i]]}

pouter6_1 <- log(pars[free_parameters6]) %>% set_names(paste0("log",names(.)))
best_fit <- myfit6 %>% as.parframe() %>% {.[2,]} %>% as.parvec() # der 2. sieht viel besser aus von der dynamik her
pouter6_1[names(best_fit)] <- best_fit

# mypred <- (g*x*p6_1)(seq(0, 48*3600, length.out = 200), pouter6_1, deriv = F)
# plotCombined(mypred, mydatalist, name%in% names(observables))

obj6_1 <- normL2(mydatalist, (g*x*p6_1))

# job6_1 <- runbg({myfit <- mstrust(objfun = obj6_1, center = pouter6_1, studyname = "methacetin", cores = 12, fits = 200, iterlim = 200, sd = 2); myfit}, machine = "knecht5", filename = "job6_1", input = global_env_without(c("mypred", "job", "myfit")))

# save(job6_1, file = "job6_1.rda")
# job6_1$check()
# load("job6_1.rda")
# myfit6_1 <- job6_1$get()$knecht5
# save(myfit6_1, file = "myfit6_1.rda")
# job6_1$purge()
load("myfit6_1.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit6_1 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit6_1 %>% as.parframe())+scale_y_log10()

mypred6_1 <- (g*x*p6_1)(mytimes*4, myfit6_1 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot6_1 <- plotCombined(mypred6_1, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot6_1)
# myplot6_1
##   index    value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1   126 2254.796      TRUE        168         -7.160945     -4.394436
## 2    55 2616.813     FALSE        200         -6.451695     -7.616642
## 3    17 2724.360      TRUE        107          5.557814      6.673239
## 4    40 3089.612      TRUE         86          4.825572      9.344937
## 5    49 3140.781     FALSE        200         -7.690501     -8.975971
## 6   152 3143.196      TRUE         57         -7.021284     -1.557577
##   logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1         -6.718847         -8.311169     -1.954417  -6.524918
## 2         -2.815882         -9.898235    -11.899718  -4.243046
## 3         -4.036551         -9.848632     -6.435657  -7.323644
## 4         -5.141421         -9.827778    -10.177033  -6.097087
## 5         -4.286810         -7.639483     -2.995771  -6.741065
## 6         -7.080733        -10.818246     -6.381171  -5.752037
##   logKpre_apap logKpki_apap logKpli_apap logKpre_apap_cys logKpki_apap_cys
## 1  0.011772621    1.1690414    0.4379879        0.9132539       -1.9017606
## 2 -0.217629833    1.1177705   -3.9170836        0.4071650        1.7266591
## 3 -0.270321088   -1.2532600   -2.0291797        0.6635614       -1.2993619
## 4 -0.296894373    2.2177932   -0.7600494       -1.0904734        2.5744728
## 5  0.005178141    1.4576398   -2.2549729       -0.7537454       -2.2248627
## 6  0.234766940    0.9271649    1.8868622       -1.9757345        0.2620245
##   logKpli_apap_cys logKpre_apap_glu logKpki_apap_glu logKpli_apap_glu
## 1       -3.2571647      -1.58858870        0.4221563       -0.4352780
## 2        0.0398897      -1.46330331        1.4233708       -0.4196139
## 3       -2.7804837       0.05884825       -0.1986642       -1.6136809
## 4        1.5410468      -3.05246767       -2.0834014       -4.2756599
## 5        2.8185400      -1.92875393        0.2224316        5.5323154
## 6       -1.5929422      -3.93035938       -0.3190701       -3.6514846
##   logKpre_apap_sul
## 1       -0.6014329
## 2       -0.4671365
## 3        0.2192896
## 4        1.7863395
## 5       -0.4661242
## 6        1.1538600
## logAPAPCYS_HLM_CL:  0.000245756553408302
##     logAPAPCYS_Km:  0.141647012033554
## logAPAPGLU_HLM_CL:  0.000776320239293284
##     logAPAPGLU_Km:  0.0123458440700416
## logAPAPSUL_HLM_CL:  0.00120793042044707
##        logKa_apap:  0.0014664396901935
##      logKpki_apap:  3.21890543991051
##  logKpki_apap_cys:  0.1493055246381
##  logKpki_apap_glu:  1.52524686894739
##      logKpli_apap:  1.54958613737682
##  logKpli_apap_cys:  0.038497394575441
##  logKpli_apap_glu:  0.647084766249551
##      logKpre_apap:  1.01184219056339
##  logKpre_apap_cys:  2.4924194943794
##  logKpre_apap_glu:  0.204213614409576
##  logKpre_apap_sul:  0.548025794974965
global_env_without <- function(reg) ls(.GlobalEnv)[!(ls(.GlobalEnv) %>% sapply(. %>% str_detect(reg) %>% any))]

Scaling factors and Error model

Introduce scaling factors 7 - not enough fits converged, but in principle not bad

load("methacetin.rda")

x <- Xs(myodemodel) # make prediction function
loadDLL(x)

# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]

free_parameters7 <- c("APAPGLU_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPSUL_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPCYS_HLM_CL",  # Vmax value
                     "APAPCYS_Km",  # Km value

                     "Ka_apap"#, #"F_apap_sul"    ,
                     # "Kpre_apap", "Kpki_apap", "Kpli_apap",
                     # "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
                     # "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
                     # "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
                     # "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
                     # "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
                     )

fixed_parameters7 <- pars[!(names(pars)%in%c(free_parameters7,names(f)[1]))] %>% names

mydatalist <- data  %>% filter(!is.na(sigma)) %>% select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")


observables7 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap", 
                 apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu", 
                 apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul", 
                 apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters7 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)

# free_parameters7 <- c(free_parameters7, scale_parameters7)
     

i <- 2
p_list <- lapply(1:nrow(conditions), function(i) {
  cond <- unlist(conditions[i,])[2:3]
  
  trafo <- as.character(pars) %>% set_names(names(pars))
  trafo[names(cond)] <- cond
  trafo[free_parameters7] <- paste0("exp(log", free_parameters7, ")")

  scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters7, x = scale_parameters7, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters7)
  scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any %>% not)] <- "1"
  
  trafo <- c(trafo, scales)
    
  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p7 <- NULL
for(i in 1:length(p_list)) { p7 <<- p7 + p_list[[i]]}

g7 <- Y(observables7, x)#, parameters = c(free_parameters7, scale_parameters7))  

obj7 <- normL2(mydatalist, (g7*x*p7))

pouter7 <- rep(0, length(getParameters(obj7))) %>% set_names(getParameters(obj7))
pouter7[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()

# job7 <- runbg({myfit <- mstrust(objfun = obj7, center = pouter7, studyname = "methacetin", cores = 12, fits = 100); myfit}, machine = "knecht5", filename = "job7", input = global_env_without(c("mypred", "job", "myfit")))

# save(job7, file = "job7.rda")
# job7$check()
# myfit7 <- job7$get()$knecht5
# save(myfit7, file = "myfit7.rda")
# job7$purge()
load("myfit7.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit7 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit7 %>% as.parframe())+scale_y_log10()

mypred7 <- (g7*x*p7)(mytimes*4, myfit7 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot7 <- plotCombined(mypred7, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot7)
# myplot7
##   index    value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1     5 2490.012      TRUE         30         -5.749303   -3.08097475
## 2    64 2626.769      TRUE         60         -4.454476   -2.31516013
## 3    67 2855.601      TRUE         87         -5.293002   -2.39553148
## 4    35 3037.778     FALSE        100         -4.055656   -2.26447288
## 5    98 3347.208     FALSE        100         -5.763790   -4.18377196
## 6    14 3572.777     FALSE        100         -2.353277   -0.00949594
##   logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1         -5.069446         -10.95827     -6.696778  -7.438818
## 2         -5.574571         -12.31058     -6.987258  -7.603783
## 3         -5.345860         -11.73765     -6.328117  -6.926927
## 4         -5.365076         -11.24999     -7.009670  -7.954771
## 5         -5.881820         -10.51708     -7.158013  -8.062787
## 6         -5.624627         -10.82433     -6.874139  -7.609710
##   log_scale_apap_Chan1997_1.4_0 log_scale_apap_glu_Chan1997_1.4_0
## 1                     0.7357347                      -0.258980454
## 2                     0.9075704                       0.003338643
## 3                     0.4662637                       0.365011207
## 4                     1.2621657                      -0.338130480
## 5                     1.0203440                      -0.713985018
## 6                     0.7625238                       0.359423877
##   log_scale_apap_sul_Chan1997_1.4_0 log_scale_apap_cys_Chan1997_1.4_0
## 1                       -0.03917457                         0.6911870
## 2                        0.21839269                         0.4441861
## 3                       -0.26185349                         1.4962657
## 4                        0.48438243                         0.9564866
## 5                        0.35501653                         0.2993361
## 6                       -0.08656187                         0.4475230
##   log_scale_apap_Chiew2010_5.6_0 log_scale_apap_glu_Chiew2010_5.6_0
## 1                     0.28907910                          1.4721450
## 2                     0.32685873                          0.8062451
## 3                     0.04167644                          0.5089211
## 4                     0.28628828                          0.1863739
## 5                    -0.25747805                         -0.5369452
## 6                     0.60731079                         -1.4552898
##   log_scale_apap_sul_Chiew2010_5.6_0 log_scale_apap_Critchley2005_1.4_0
## 1                         -0.1620631                         -0.3123042
## 2                          0.8175662                          0.1778291
## 3                         -2.1840586                         -0.8736844
## 4                          0.8380141                          0.2651291
## 5                          0.2400412                         -1.6486526
## 6                          0.6459171                          0.1341244
##   log_scale_apap_glu_Critchley2005_1.4_0
## 1                              0.6340540
## 2                             -0.8673394
## 3                             -0.5829681
## 4                              0.1720526
## 5                             -0.6356093
## 6                             -1.1555520
##   log_scale_apap_sul_Critchley2005_1.4_0
## 1                            -0.32952950
## 2                            -0.58790653
## 3                            -0.75363120
## 4                            -0.42899409
## 5                             0.36572246
## 6                             0.06592612
##   log_scale_apap_cys_Critchley2005_1.4_0 log_scale_apap_Rawlins1977_0_1
## 1                              0.4520727                     0.12354755
## 2                              1.9508658                     0.19603601
## 3                              1.1278751                     0.01086081
## 4                              0.8794284                     0.34127665
## 5                             -0.4761007                     0.23294843
## 6                              0.3483191                     0.15601077
##   log_scale_apap_Rawlins1977_1_0 log_scale_apap_Rawlins1977_2_0
## 1                     0.27999826                     0.23205478
## 2                     0.40194965                     0.40054980
## 3                     0.07005591                     0.03746354
## 4                     0.69124502                     0.72483500
## 5                     0.55520067                     0.47151149
## 6                     0.28847757                     0.35586546
##   log_scale_apap_Rawlins1977_0.5_0
## 1                      -0.32097913
## 2                      -0.03674230
## 3                      -0.03173343
## 4                      -0.78815698
## 5                       0.17669530
## 6                       0.45105346
##                      logAPAPCYS_HLM_CL:  1.74133804776165e-05
##                          logAPAPCYS_Km:  0.00123488444599134
##                      logAPAPGLU_HLM_CL:  0.0031850000915229
##                          logAPAPGLU_Km:  0.0459144794614278
##                      logAPAPSUL_HLM_CL:  0.00628590132812984
##                             logKa_apap:  0.000587979808605786
##          log_scale_apap_Chan1997_1.4_0:  2.08701470180148
##         log_scale_apap_Chiew2010_5.6_0:  1.3351973371447
##     log_scale_apap_Critchley2005_1.4_0:  0.731758900597009
##      log_scale_apap_cys_Chan1997_1.4_0:  1.99608345214162
## log_scale_apap_cys_Critchley2005_1.4_0:  1.57156625159931
##      log_scale_apap_glu_Chan1997_1.4_0:  0.771838109420778
##     log_scale_apap_glu_Chiew2010_5.6_0:  4.35857412292258
## log_scale_apap_glu_Critchley2005_1.4_0:  1.8852379216948
##         log_scale_apap_Rawlins1977_0_1:  1.13150380172748
##       log_scale_apap_Rawlins1977_0.5_0:  0.725438390393933
##         log_scale_apap_Rawlins1977_1_0:  1.32312750633397
##         log_scale_apap_Rawlins1977_2_0:  1.26118881349357
##      log_scale_apap_sul_Chan1997_1.4_0:  0.96158283149328
##     log_scale_apap_sul_Chiew2010_5.6_0:  0.850387509055339
## log_scale_apap_sul_Critchley2005_1.4_0:  0.719262064522545

Error model

8 Introduce error model and scaling factors in the dynamic model - very weird results.

Very bad results… Mistake: Don’t remove the rows with sigma=NA! Stupid…

load("methacetin.rda")

x <- Xs(myodemodel) # make prediction function
loadDLL(x)

# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]

free_parameters8 <- c("APAPGLU_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPSUL_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPCYS_HLM_CL",  # Vmax value
                     "APAPCYS_Km",  # Km value

                     "Ka_apap"#, #"F_apap_sul"    ,
                     # "Kpre_apap", "Kpki_apap", "Kpli_apap",
                     # "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
                     # "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
                     # "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
                     # "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
                     # "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
                     )

fixed_parameters8 <- pars[!(names(pars)%in%c(free_parameters8,names(f)[1]))] %>% names

mydatalist <- data  %>% filter(!is.na(sigma)) %>% select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")


observables8 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap", 
                 apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu", 
                 apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul", 
                 apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters8 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)

# free_parameters8 <- c(free_parameters8, scale_parameters8)

error_model8 <- c(apap = "srel_apap*apap^2 +s0_apap", 
                 apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu", 
                 apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul", 
                 apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
error_parameters8 <- setdiff(getSymbols(error_model8), names(error_model8)) %>% set_names(.,.)

i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
  cond <- unlist(conditions[i,])[2:3]
  
  trafo <- as.character(pars) %>% set_names(names(pars))
  trafo[names(cond)] <- cond
  trafo[free_parameters8] <- paste0("exp(log", free_parameters8, ")")

  scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters8, x = scale_parameters8, y = .)} 
  scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  
  errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters8, x = error_parameters8, y = .)}
  errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  
  trafo <- c(trafo, scales, errors)
    
  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p8 <- NULL
for(i in 1:length(p_list)) { p8 <<- p8 + p_list[[i]]}

g8 <- Y(observables8, x)#, parameters = c(free_parameters8, scale_parameters8))  

err8 <- Y(error_model8, g8)

obj8 <- normL2(mydatalist, (g8*x*p8), errmodel = err8)

pouter8 <- rep(0, length(getParameters(obj8))) %>% set_names(getParameters(obj8))
pouter8[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()

# obj8(pouter8)

# job8 <- runbg({myfit <- mstrust(objfun = obj8, center = pouter8, studyname = "methacetin", cores = 12, fits = 100); myfit}, machine = "knecht5", filename = "job8", input = global_env_without(c("mypred", "job", "myfit")))

# save(job8, file = "job8.rda")

# job8$purge()
# myfit8 <- job8$get()$knecht5
# save(myfit8, file = "myfit8.rda")
# job8$purge()
load("myfit8.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit8 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit8 %>% as.parframe())#+scale_y_log10()

mypred8 <- (g8*x*p8)(mytimes*4, myfit8 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot8 <- plotCombined(mypred8, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot8)
# myplot8
##   index     value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1    46 -486.5434      TRUE         68         0.8448040     3.0690138
## 2    95 -477.9974      TRUE         86         2.3221006     4.5447823
## 3     3 -458.5112      TRUE         85        -2.8788000    -0.5052364
## 4    91 -451.1060      TRUE         84         2.0359641     4.2579075
## 5    50 -437.6416     FALSE        100         2.1876617     4.4051112
## 6    44 -426.7063      TRUE         82        -0.9089178     1.3108703
##   logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1        -10.788256         -13.56917     -3.126300  -6.295179
## 2        -12.281700         -11.21214     -1.128100  -6.280291
## 3         -6.651499         -14.82540     -3.552027  -6.325155
## 4        -12.473814         -13.54999     -2.420064  -6.289021
## 5        -11.711911         -13.73429     -2.595443  -2.278541
## 6        -10.442107         -14.43848     -3.406845  -2.307167
##   log_scale_apap_Chan1997_1.4_0 log_scale_apap_glu_Chan1997_1.4_0
## 1                   -0.34147420                         0.3524540
## 2                    0.87395974                        -0.2202193
## 3                    0.28419894                        -1.2117128
## 4                   -0.23704451                        -0.5141464
## 5                   -0.00256986                         0.7640506
## 6                   -1.39510256                        -0.8312082
##   log_scale_apap_sul_Chan1997_1.4_0 log_scale_apap_cys_Chan1997_1.4_0
## 1                       -0.09054225                         0.4909674
## 2                       -2.05176756                         0.1713548
## 3                        1.59152206                         0.9567482
## 4                       -0.64296093                         0.6822799
## 5                       -0.28381384                         1.2033634
## 6                       -0.42447087                        -0.9383802
##   log_srel_apap_Chan1997_1.4_0 log_s0_apap_Chan1997_1.4_0
## 1                    0.2569589                 -1.9014794
## 2                   -0.6194964                  1.4224062
## 3                    0.4884257                  0.4620627
## 4                    0.4431686                  0.3437440
## 5                   -1.0774869                 -0.7111688
## 6                   -0.1059288                 -0.4287323
##   log_srel_apap_glu_Chan1997_1.4_0 log_s0_apap_glu_Chan1997_1.4_0
## 1                       -0.9569274                     -0.1547545
## 2                        0.6830946                     -0.9358321
## 3                        0.1252113                      1.4079458
## 4                       -1.0505252                     -1.1969937
## 5                        1.3007331                     -1.5656455
## 6                        2.1081213                     -1.7762023
##   log_srel_apap_sul_Chan1997_1.4_0 log_s0_apap_sul_Chan1997_1.4_0
## 1                       -2.5583855                     -1.1857697
## 2                        0.2056214                     -1.4048882
## 3                       -0.5206815                     -3.0952762
## 4                        1.7018681                     -2.4248608
## 5                       -1.9290153                     -0.8441127
## 6                        0.6754411                      0.6658790
##   log_srel_apap_cys_Chan1997_1.4_0 log_s0_apap_cys_Chan1997_1.4_0
## 1                      -0.80413329                    -1.76384101
## 2                      -1.03467092                    -2.63912882
## 3                       0.12948514                    -0.06831229
## 4                      -1.23494485                    -0.35465510
## 5                      -1.45640160                    -2.21684171
## 6                      -0.08965077                    -2.18967728
##   log_scale_apap_Chiew2010_5.6_0 log_scale_apap_glu_Chiew2010_5.6_0
## 1                    -0.73225919                         0.29875836
## 2                     0.07746059                        -0.75638815
## 3                    -0.41557105                         0.66860572
## 4                    -0.34446349                        -0.09935669
## 5                    -0.06790696                         0.17262465
## 6                    -0.36695303                         0.98236576
##   log_scale_apap_sul_Chiew2010_5.6_0 log_srel_apap_Chiew2010_5.6_0
## 1                         -0.2751353                     0.2001658
## 2                         -0.8853731                    -1.1763805
## 3                         -1.4579817                     1.4459862
## 4                          1.9013839                    -1.1825723
## 5                          0.8957957                     0.4329237
## 6                         -0.5285034                    -1.6273736
##   log_s0_apap_Chiew2010_5.6_0 log_srel_apap_glu_Chiew2010_5.6_0
## 1                  -1.3083801                         0.2537679
## 2                  -0.9910735                         0.4781154
## 3                  -1.3495770                        -0.2456441
## 4                  -0.4922068                        -0.4925589
## 5                   1.2760837                         1.4551869
## 6                  -1.4419611                        -2.2427514
##   log_s0_apap_glu_Chiew2010_5.6_0 log_srel_apap_sul_Chiew2010_5.6_0
## 1                      -0.9917074                        0.30899346
## 2                      -0.1322426                        0.20086298
## 3                      -2.4500674                        1.13602514
## 4                      -2.1513778                        0.97349583
## 5                       0.3285409                        0.09686059
## 6                      -1.1718316                       -0.32388955
##   log_s0_apap_sul_Chiew2010_5.6_0 log_scale_apap_Critchley2005_1.4_0
## 1                      -0.2317920                         0.75661973
## 2                      -0.1080672                         0.40330247
## 3                      -1.2538891                        -1.84876130
## 4                       0.5970107                         0.09984706
## 5                       0.3377952                        -1.08319140
## 6                      -0.9579787                         0.15268520
##   log_scale_apap_glu_Critchley2005_1.4_0
## 1                             -1.6075093
## 2                             -0.4203937
## 3                             -0.4089714
## 4                              0.6322839
## 5                              0.5396460
## 6                              0.6909803
##   log_scale_apap_sul_Critchley2005_1.4_0
## 1                              1.1200586
## 2                              1.3772116
## 3                             -0.8219418
## 4                             -0.1005233
## 5                             -1.1104974
## 6                             -0.4086092
##   log_scale_apap_cys_Critchley2005_1.4_0 log_srel_apap_Critchley2005_1.4_0
## 1                              0.7914627                         0.5229484
## 2                             -0.2206032                         1.7093085
## 3                             -0.7028494                        -0.4162756
## 4                             -0.3442296                        -1.5224635
## 5                             -0.7764242                        -1.5065604
## 6                             -1.2219040                        -0.4096500
##   log_s0_apap_Critchley2005_1.4_0 log_srel_apap_glu_Critchley2005_1.4_0
## 1                       2.8885995                            -0.6958732
## 2                       0.2769623                             1.6113013
## 3                       1.9668600                             0.8816571
## 4                       0.4696681                             0.6190646
## 5                      -0.1120432                             0.9772918
## 6                       2.0007139                            -1.2027400
##   log_s0_apap_glu_Critchley2005_1.4_0
## 1                        -1.462525446
## 2                        -0.876828265
## 3                        -2.072184801
## 4                        -0.374247260
## 5                        -0.093214581
## 6                         0.005710159
##   log_srel_apap_sul_Critchley2005_1.4_0
## 1                            0.09474047
## 2                            1.05831194
## 3                            1.37006018
## 4                            0.49068959
## 5                           -0.15389290
## 6                           -0.21423115
##   log_s0_apap_sul_Critchley2005_1.4_0
## 1                          -1.5436022
## 2                          -1.9787239
## 3                           1.1611746
## 4                          -0.5722336
## 5                          -0.7604901
## 6                          -0.5263634
##   log_srel_apap_cys_Critchley2005_1.4_0
## 1                             0.6346692
## 2                            -0.9788842
## 3                            -1.7831523
## 4                             0.6156580
## 5                            -1.3146436
## 6                             0.4575737
##   log_s0_apap_cys_Critchley2005_1.4_0 log_scale_apap_Rawlins1977_0_1
## 1                           0.4977260                    0.021621889
## 2                          -0.6198171                    0.021845011
## 3                          -1.8375393                    0.008699087
## 4                          -1.2241609                    0.022133087
## 5                          -1.1091981                    0.026012429
## 6                          -0.1110733                    0.024257471
##   log_srel_apap_Rawlins1977_0_1 log_s0_apap_Rawlins1977_0_1
## 1                      2.454559                   -8.977730
## 2                      2.454535                   -8.977814
## 3                      2.483257                   -8.996088
## 4                      2.453265                   -8.974551
## 5                      2.437694                   -8.941497
## 6                      2.442176                   -8.947284
##   log_scale_apap_Rawlins1977_1_0 log_srel_apap_Rawlins1977_1_0
## 1                      0.1678273                      2.582696
## 2                      0.1692129                      2.581302
## 3                      0.1388991                      2.595842
## 4                      0.1692214                      2.580949
## 5                      0.2307545                      2.589881
## 6                      0.2203842                      2.591367
##   log_s0_apap_Rawlins1977_1_0 log_scale_apap_Rawlins1977_2_0
## 1                   -14.66563                      0.4137086
## 2                   -16.05141                      0.4157914
## 3                   -15.53460                      0.3500338
## 4                   -15.71205                      0.4157741
## 5                   -16.02509                      0.4775457
## 6                   -15.33039                      0.4570333
##   log_srel_apap_Rawlins1977_2_0 log_s0_apap_Rawlins1977_2_0
## 1                      2.341414                   -15.33500
## 2                      2.341451                   -16.63546
## 3                      2.347062                   -16.04883
## 4                      2.341416                   -16.30266
## 5                      2.346485                   -16.51624
## 6                      2.347016                   -15.82922
##   log_scale_apap_Rawlins1977_0.5_0 log_srel_apap_Rawlins1977_0.5_0
## 1                       -1.5900253                       1.9345084
## 2                       -0.2137423                       0.2236760
## 3                        0.8627395                       0.6101951
## 4                        1.1149818                      -0.5099961
## 5                       -0.4089437                       1.0647446
## 6                       -0.1155511                       0.7124745
##   log_s0_apap_Rawlins1977_0.5_0
## 1                    -2.5199219
## 2                    -0.2345078
## 3                    -1.2819850
## 4                     1.5677611
## 5                    -0.7884726
## 6                    -0.3613277
##                      logAPAPCYS_HLM_CL:  1.2793391679861e-06
##                          logAPAPCYS_Km:  0.0438798488903149
##                      logAPAPGLU_HLM_CL:  2.32752164063581
##                          logAPAPGLU_Km:  21.520668349057
##                      logAPAPSUL_HLM_CL:  2.06404786189464e-05
##                             logKa_apap:  0.00184517865281865
##             log_s0_apap_Chan1997_1.4_0:  0.149347509305303
##            log_s0_apap_Chiew2010_5.6_0:  0.270257502535883
##        log_s0_apap_Critchley2005_1.4_0:  17.9681281235639
##         log_s0_apap_cys_Chan1997_1.4_0:  0.171385306079452
##    log_s0_apap_cys_Critchley2005_1.4_0:  1.64497626822607
##         log_s0_apap_glu_Chan1997_1.4_0:  0.856625476079955
##        log_s0_apap_glu_Chiew2010_5.6_0:  0.370942798641939
##    log_s0_apap_glu_Critchley2005_1.4_0:  0.231650514577204
##            log_s0_apap_Rawlins1977_0_1:  0.000126188936947135
##          log_s0_apap_Rawlins1977_0.5_0:  0.0804658932475238
##            log_s0_apap_Rawlins1977_1_0:  4.27362459933834e-07
##            log_s0_apap_Rawlins1977_2_0:  2.18822833245219e-07
##         log_s0_apap_sul_Chan1997_1.4_0:  0.30551094502221
##        log_s0_apap_sul_Chiew2010_5.6_0:  0.793111039544634
##    log_s0_apap_sul_Critchley2005_1.4_0:  0.21361025170362
##          log_scale_apap_Chan1997_1.4_0:  0.710721805251837
##         log_scale_apap_Chiew2010_5.6_0:  0.480821497303708
##     log_scale_apap_Critchley2005_1.4_0:  2.13106047495752
##      log_scale_apap_cys_Chan1997_1.4_0:  1.63389603961365
## log_scale_apap_cys_Critchley2005_1.4_0:  2.20662164565056
##      log_scale_apap_glu_Chan1997_1.4_0:  1.42255421536611
##     log_scale_apap_glu_Chiew2010_5.6_0:  1.34818380327456
## log_scale_apap_glu_Critchley2005_1.4_0:  0.200386104365974
##         log_scale_apap_Rawlins1977_0_1:  1.02185733567333
##       log_scale_apap_Rawlins1977_0.5_0:  0.203920455093923
##         log_scale_apap_Rawlins1977_1_0:  1.18273235599321
##         log_scale_apap_Rawlins1977_2_0:  1.51241638200544
##      log_scale_apap_sul_Chan1997_1.4_0:  0.913435744877439
##     log_scale_apap_sul_Chiew2010_5.6_0:  0.759469382839319
## log_scale_apap_sul_Critchley2005_1.4_0:  3.06503365611115
##           log_srel_apap_Chan1997_1.4_0:  1.29299192947592
##          log_srel_apap_Chiew2010_5.6_0:  1.22160522476788
##      log_srel_apap_Critchley2005_1.4_0:  1.68699428030717
##       log_srel_apap_cys_Chan1997_1.4_0:  0.447475592070443
##  log_srel_apap_cys_Critchley2005_1.4_0:  1.88639794692755
##       log_srel_apap_glu_Chan1997_1.4_0:  0.384071160069748
##      log_srel_apap_glu_Chiew2010_5.6_0:  1.28887257124136
##  log_srel_apap_glu_Critchley2005_1.4_0:  0.498638857365226
##          log_srel_apap_Rawlins1977_0_1:  11.6412981287543
##        log_srel_apap_Rawlins1977_0.5_0:  6.92064127305638
##          log_srel_apap_Rawlins1977_1_0:  13.2327664425859
##          log_srel_apap_Rawlins1977_2_0:  10.3959247377831
##       log_srel_apap_sul_Chan1997_1.4_0:  0.0774296530505137
##      log_srel_apap_sul_Chiew2010_5.6_0:  1.36205346050779
##  log_srel_apap_sul_Critchley2005_1.4_0:  1.09937349942295

9 Introduce error model and scaling factors in the dynamic model - not so bad, but not good either

  • some scaling factors look as if they could help
  • dynamics is described quite ok
  • in apap_sul, model is still too slow
  • in apap_glu, model is a bit too slow
  • chiew 2010, apap_glu is better than in try 10.
load("methacetin.rda")

x <- Xs(myodemodel) # make prediction function
loadDLL(x)

# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]

free_parameters9 <- c("APAPGLU_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPSUL_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPCYS_HLM_CL",  # Vmax value
                     "APAPCYS_Km",  # Km value

                     "Ka_apap"#, #"F_apap_sul"    ,
                     # "Kpre_apap", "Kpki_apap", "Kpli_apap",
                     # "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
                     # "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
                     # "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
                     # "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
                     # "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
                     )

fixed_parameters9 <- pars[!(names(pars)%in%c(free_parameters9,names(f)[1]))] %>% names

mydatalist <- data  %>% select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")


observables9 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap", 
                 apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu", 
                 apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul", 
                 apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters9 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)

# free_parameters9 <- c(free_parameters9, scale_parameters9)

error_model9 <- c(apap = "srel_apap*apap^2 +s0_apap", 
                 apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu", 
                 apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul", 
                 apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
error_parameters9 <- setdiff(getSymbols(error_model9), names(error_model9)) %>% set_names(.,.)

i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
  cond <- unlist(conditions[i,])[2:3]
  
  trafo <- as.character(pars) %>% set_names(names(pars))
  trafo[names(cond)] <- cond
  trafo[free_parameters9] <- paste0("exp(log", free_parameters9, ")")

  scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters9, x = scale_parameters9, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters9)
  scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  
  errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters9, x = error_parameters9, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters9)
  errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  
  trafo <- c(trafo, scales, errors)
    
  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p9 <- NULL
for(i in 1:length(p_list)) { p9 <<- p9 + p_list[[i]]}

g9 <- Y(observables9, x)#, parameters = c(free_parameters9, scale_parameters9))  

err9 <- Y(error_model9, g9)

obj9 <- normL2(mydatalist, (g9*x*p9), errmodel = err9)

pouter9 <- rep(0, length(getParameters(obj9))) %>% set_names(getParameters(obj9))
pouter9[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()

# obj9(pouter9)

# job9 <- runbg({myfit <- mstrust(objfun = obj9, center = pouter9, studyname = "methacetin", cores = 16, fits = 200); myfit}, machine = "ruprecht2", filename = "job9", input = global_env_without(c("mypred", "job", "myfit")))

# save(job9, file = "job9.rda")

# job9$check()
# 
# load("job9.rda")
# myfit9 <- job9$get()$ruprecht2
# save(myfit9, file = "myfit9.rda")
# job9$purge()
load("myfit9.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit9 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit9 %>% as.parframe())#+scale_y_log10()

mypred9 <- (g9*x*p9)(mytimes*4, myfit9 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot9 <- plotCombined(mypred9, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot9)
# myplot9
# obj9(myfit9 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.})
##   index     value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1     1 -2477.716     FALSE        100         -5.099331      7.999288
## 2    35 -2477.716     FALSE        100         -5.032580      7.301881
## 3    49 -2477.715     FALSE        100         -5.087414      7.628175
## 4   179 -2477.715     FALSE        100         -5.087632      7.116519
## 5    64 -2477.715      TRUE         99         -5.073408      7.132834
## 6   127 -2477.715     FALSE        100         -5.094321      6.967860
##   logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1         -4.797299         -17.83559     -6.250228  -7.073747
## 2         -4.797371         -17.08140     -6.250285  -7.073740
## 3         -4.797351         -17.45008     -6.250078  -7.073678
## 4         -4.797346         -16.93769     -6.250345  -7.073796
## 5         -4.797347         -16.94634     -6.250387  -7.073790
## 6         -4.797427         -16.79246     -6.250201  -7.073706
##   log_scale_apap_Chan1997_14_0 log_scale_apap_glu_Chan1997_14_0
## 1                    0.3201883                        10.240255
## 2                    0.3201823                         9.476110
## 3                    0.3201498                         9.857194
## 4                    0.3202139                         9.345829
## 5                    0.3202098                         9.347900
## 6                    0.3201586                         9.203812
##   log_scale_apap_sul_Chan1997_14_0 log_scale_apap_cys_Chan1997_14_0
## 1                       -0.5703412                         7.541992
## 2                       -0.5702552                         6.787785
## 3                       -0.5703223                         7.156501
## 4                       -0.5702481                         6.644072
## 5                       -0.5702448                         6.652709
## 6                       -0.5702149                         6.498853
##   log_srel_apap_Chan1997_14_0 log_s0_apap_Chan1997_14_0
## 1                   -1.066589                 -8.553130
## 2                   -1.066910                 -8.552872
## 3                   -1.064453                 -8.553918
## 4                   -1.069107                 -8.552475
## 5                   -1.070049                 -8.552360
## 6                   -1.066075                 -8.552996
##   log_srel_apap_glu_Chan1997_14_0 log_s0_apap_glu_Chan1997_14_0
## 1                        3.937521                     -7.368639
## 2                        3.937853                     -7.368713
## 3                        3.937651                     -7.368663
## 4                        3.937780                     -7.368695
## 5                        3.937766                     -7.368738
## 6                        3.937715                     -7.368682
##   log_srel_apap_sul_Chan1997_14_0 log_s0_apap_sul_Chan1997_14_0
## 1                        4.178282                     -9.337702
## 2                        4.178646                     -9.337195
## 3                        4.178415                     -9.337887
## 4                        4.178152                     -9.337426
## 5                        4.178408                     -9.337071
## 6                        4.178751                     -9.337377
##   log_srel_apap_cys_Chan1997_14_0 log_s0_apap_cys_Chan1997_14_0
## 1                       -4.014309                     -9.955108
## 2                       -3.320616                     -9.955241
## 3                       -3.635824                     -9.955119
## 4                       -3.161390                     -9.955246
## 5                       -3.162397                     -9.955206
## 6                       -2.983036                     -9.955250
##   log_scale_apap_Chiew2010_56_0 log_scale_apap_glu_Chiew2010_56_0
## 1                     0.1542391                         10.561511
## 2                     0.1542455                          9.797452
## 3                     0.1541823                         10.178496
## 4                     0.1542872                          9.667135
## 5                     0.1542769                          9.669024
## 6                     0.1542120                          9.525162
##   log_scale_apap_sul_Chiew2010_56_0 log_srel_apap_Chiew2010_56_0
## 1                        -0.7426937                    0.3778347
## 2                        -0.7426260                    0.3779957
## 3                        -0.7426694                    0.3771786
## 4                        -0.7426272                    0.3780723
## 5                        -0.7426270                    0.3777571
## 6                        -0.7425874                    0.3779534
##   log_s0_apap_Chiew2010_56_0 log_srel_apap_glu_Chiew2010_56_0
## 1                  -7.069136                         2.369965
## 2                  -7.067144                         2.370166
## 3                  -7.067803                         2.370305
## 4                  -7.066808                         2.371076
## 5                  -7.066785                         2.371267
## 6                  -7.067703                         2.370746
##   log_s0_apap_glu_Chiew2010_56_0 log_srel_apap_sul_Chiew2010_56_0
## 1                      -5.863964                         1.413941
## 2                      -5.863894                         1.413643
## 3                      -5.864205                         1.414305
## 4                      -5.864437                         1.413383
## 5                      -5.865029                         1.413175
## 6                      -5.864204                         1.413946
##   log_s0_apap_sul_Chiew2010_56_0 log_scale_apap_Critchley2005_14_0
## 1                      -6.186949                         0.2345061
## 2                      -6.187438                         0.2344554
## 3                      -6.187013                         0.2343445
## 4                      -6.187419                         0.2344923
## 5                      -6.187215                         0.2345560
## 6                      -6.187368                         0.2344314
##   log_scale_apap_glu_Critchley2005_14_0
## 1                             10.314313
## 2                              9.550148
## 3                              9.931238
## 4                              9.419870
## 5                              9.421957
## 6                              9.277840
##   log_scale_apap_sul_Critchley2005_14_0
## 1                             -1.148680
## 2                             -1.148614
## 3                             -1.148658
## 4                             -1.148615
## 5                             -1.148612
## 6                             -1.148576
##   log_scale_apap_cys_Critchley2005_14_0 log_srel_apap_Critchley2005_14_0
## 1                              7.232789                         2.288037
## 2                              6.478583                         2.284297
## 3                              6.847288                         2.279003
## 4                              6.334869                         2.284317
## 5                              6.343503                         2.285165
## 6                              6.189646                         2.284663
##   log_s0_apap_Critchley2005_14_0 log_srel_apap_glu_Critchley2005_14_0
## 1                      -8.080480                            -5.922499
## 2                      -8.074298                            -5.172423
## 3                      -8.068500                            -5.541084
## 4                      -8.073987                            -5.031696
## 5                      -8.075026                            -5.034967
## 6                      -8.075041                            -4.887071
##   log_s0_apap_glu_Critchley2005_14_0 log_srel_apap_sul_Critchley2005_14_0
## 1                          -6.304482                             4.463234
## 2                          -6.304953                             4.464450
## 3                          -6.304470                             4.464503
## 4                          -6.304871                             4.464421
## 5                          -6.304808                             4.464391
## 6                          -6.304898                             4.464320
##   log_s0_apap_sul_Critchley2005_14_0 log_srel_apap_cys_Critchley2005_14_0
## 1                          -7.946701                            -2.811800
## 2                          -7.945526                            -2.085599
## 3                          -7.945487                            -2.428260
## 4                          -7.945550                            -1.991045
## 5                          -7.945590                            -1.941678
## 6                          -7.945679                            -1.917429
##   log_s0_apap_cys_Critchley2005_14_0 log_scale_apap_Rawlins1977_0_1
## 1                          -10.12659                     -0.1015749
## 2                          -10.12657                     -0.1015765
## 3                          -10.12661                     -0.1016133
## 4                          -10.12658                     -0.1015472
## 5                          -10.12661                     -0.1015634
## 6                          -10.12657                     -0.1015950
##   log_srel_apap_Rawlins1977_0_1 log_s0_apap_Rawlins1977_0_1
## 1                      2.881313                   -10.45464
## 2                      2.880629                   -10.44504
## 3                      2.880721                   -10.44162
## 4                      2.880528                   -10.44397
## 5                      2.879162                   -10.41924
## 6                      2.880092                   -10.43567
##   log_scale_apap_Rawlins1977_1_0 log_srel_apap_Rawlins1977_1_0
## 1                    -0.03506085                      3.061212
## 2                    -0.03507190                      3.062111
## 3                    -0.03511765                      3.062435
## 4                    -0.03503499                      3.062118
## 5                    -0.03503313                      3.062624
## 6                    -0.03510720                      3.062807
##   log_s0_apap_Rawlins1977_1_0 log_scale_apap_Rawlins1977_2_0
## 1                   -8.476342                     0.03209758
## 2                   -8.477932                     0.03208771
## 3                   -8.478216                     0.03205411
## 4                   -8.478040                     0.03212816
## 5                   -8.478879                     0.03214222
## 6                   -8.479337                     0.03205324
##   log_srel_apap_Rawlins1977_2_0 log_s0_apap_Rawlins1977_2_0
## 1                      2.456148                   -8.369891
## 2                      2.457079                   -8.371340
## 3                      2.457015                   -8.370964
## 4                      2.457062                   -8.371123
## 5                      2.457133                   -8.372795
## 6                      2.456935                   -8.370860
##   log_scale_apap_Rawlins1977_05_0 log_srel_apap_Rawlins1977_05_0
## 1                      -0.2098564                      -7.088126
## 2                      -0.2098571                      -6.338162
## 3                      -0.2098887                      -6.706408
## 4                      -0.2098285                      -6.195749
## 5                      -0.2098314                      -6.201175
## 6                      -0.2098756                      -6.053133
##   log_s0_apap_Rawlins1977_05_0
## 1                    -8.283322
## 2                    -8.284305
## 3                    -8.283902
## 4                    -8.284432
## 5                    -8.284290
## 6                    -8.284218
##                     logAPAPCYS_HLM_CL:  1.795146994404e-08
##                         logAPAPCYS_Km:  0.00193001355192211
##                     logAPAPGLU_HLM_CL:  0.00610082435223236
##                         logAPAPGLU_Km:  2978.83487681056
##                     logAPAPSUL_HLM_CL:  0.00825200823079057
##                            logKa_apap:  0.000847053225794413
##             log_s0_apap_Chan1997_14_0:  0.00019294018495771
##            log_s0_apap_Chiew2010_56_0:  0.00085096838208337
##        log_s0_apap_Critchley2005_14_0:  0.000309522401143699
##         log_s0_apap_cys_Chan1997_14_0:  4.74844786235421e-05
##    log_s0_apap_cys_Critchley2005_14_0:  4.00018130602822e-05
##         log_s0_apap_glu_Chan1997_14_0:  0.000630726272006818
##        log_s0_apap_glu_Chiew2010_56_0:  0.00283996453762311
##    log_s0_apap_glu_Critchley2005_14_0:  0.00182809286937598
##           log_s0_apap_Rawlins1977_0_1:  2.8814215684246e-05
##          log_s0_apap_Rawlins1977_05_0:  0.000252696390036355
##           log_s0_apap_Rawlins1977_1_0:  0.000208339338154456
##           log_s0_apap_Rawlins1977_2_0:  0.000231740837081231
##         log_s0_apap_sul_Chan1997_14_0:  8.80414871966006e-05
##        log_s0_apap_sul_Chiew2010_56_0:  0.00205609044641647
##    log_s0_apap_sul_Critchley2005_14_0:  0.000353827456947821
##          log_scale_apap_Chan1997_14_0:  1.37738706359055
##         log_scale_apap_Chiew2010_56_0:  1.16676983632714
##     log_scale_apap_Critchley2005_14_0:  1.26428417568878
##      log_scale_apap_cys_Chan1997_14_0:  1885.5826917627
## log_scale_apap_cys_Critchley2005_14_0:  1384.0773914102
##      log_scale_apap_glu_Chan1997_14_0:  28008.2659198081
##     log_scale_apap_glu_Chiew2010_56_0:  38619.4394146001
## log_scale_apap_glu_Critchley2005_14_0:  30161.2433788063
##        log_scale_apap_Rawlins1977_0_1:  0.903413505390761
##       log_scale_apap_Rawlins1977_05_0:  0.81070064561698
##        log_scale_apap_Rawlins1977_1_0:  0.965546657143368
##        log_scale_apap_Rawlins1977_2_0:  1.03261825838704
##      log_scale_apap_sul_Chan1997_14_0:  0.565332533789827
##     log_scale_apap_sul_Chiew2010_56_0:  0.475830435424811
## log_scale_apap_sul_Critchley2005_14_0:  0.317054858672663
##           log_srel_apap_Chan1997_14_0:  0.344180473770218
##          log_srel_apap_Chiew2010_56_0:  1.45912166820889
##      log_srel_apap_Critchley2005_14_0:  9.85557305660906
##       log_srel_apap_cys_Chan1997_14_0:  0.0180554247494083
##  log_srel_apap_cys_Critchley2005_14_0:  0.0600966934703123
##       log_srel_apap_glu_Chan1997_14_0:  51.2912931730463
##      log_srel_apap_glu_Chiew2010_56_0:  10.6970164613051
##  log_srel_apap_glu_Critchley2005_14_0:  0.00267849706185607
##         log_srel_apap_Rawlins1977_0_1:  17.8376851856982
##        log_srel_apap_Rawlins1977_05_0:  0.000834960350013315
##         log_srel_apap_Rawlins1977_1_0:  21.3534269462142
##         log_srel_apap_Rawlins1977_2_0:  11.6598121177031
##       log_srel_apap_sul_Chan1997_14_0:  65.2536831006089
##      log_srel_apap_sul_Chiew2010_56_0:  4.11213078704973
##  log_srel_apap_sul_Critchley2005_14_0:  86.7676957811533

10 Introduce error model only - quite nice, problems in chiew2100 and in apap_sul

  • some scaling factors look as if they could help
  • dynamics is described quite ok
  • in apap_sul, model is still too slow
  • in apap_glu, model is a bit too slow
  • chiew 2010, apap_sul is better than in tr9
load("methacetin.rda")

x <- Xs(myodemodel) # make prediction function
loadDLL(x)

# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]

free_parameters10 <- c("APAPGLU_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPSUL_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPCYS_HLM_CL",  # Vmax value
                     "APAPCYS_Km",  # Km value

                     "Ka_apap"#, #"F_apap_sul"    ,
                     # "Kpre_apap", "Kpki_apap", "Kpli_apap",
                     # "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
                     # "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
                     # "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
                     # "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
                     # "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
                     )

fixed_parameters10 <- pars[!(names(pars)%in%c(free_parameters10,names(f)[1]))] %>% names

mydatalist <- data %>% select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")


observables10 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
                 apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
                 apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
                 apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters10 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)

# free_parameters10 <- c(free_parameters10, scale_parameters10)

error_model10 <- c(apap = "srel_apap*apap^2 +s0_apap", 
                 apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu", 
                 apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul", 
                 apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
error_parameters10 <- setdiff(getSymbols(error_model10), names(error_model10)) %>% set_names(.,.)

p_list <- lapply(1:nrow(conditions), function(i) {
  cond <- unlist(conditions[i,])[2:3]
  
  trafo <- as.character(pars) %>% set_names(names(pars))
  trafo[names(cond)] <- cond
  trafo[free_parameters10] <- paste0("exp(log", free_parameters10, ")")

  scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters10, x = scale_parameters10, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters10)
  scales[1:length(scales)] <- "1"
  
  errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters10, x = error_parameters10, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters10)
  errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  
  trafo <- c(trafo, scales, errors)
    
  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p10 <- NULL
for(i in 1:length(p_list)) { p10 <<- p10 + p_list[[i]]}

g10 <- Y(observables10, x)#, parameters = c(free_parameters10, scale_parameters10))  

err10 <- Y(error_model10, g10)

obj10 <- normL2(mydatalist, (g10*x*p10), errmodel = err10)


pouter10 <- rep(0, length(getParameters(obj10))) %>% set_names(getParameters(obj10))
pouter10[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()

# job10 <- runbg({myfit <- mstrust(objfun = obj10, center = pouter10, studyname = "methacetin", cores = 16, fits = 200); myfit}, machine = "ruprecht1", filename = "job10", input = global_env_without(c("mypred", "job", "myfit")))

# save(job10, file = "job10.rda")

# job10$check()
# load("job10.rda")
# myfit10 <- job10$get()$ruprecht1
# save(myfit10, file = "myfit10.rda")
# job10$purge()
load("myfit10.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit10 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit10 %>% as.parframe())#+scale_y_log10()

mypred10 <- (g10*x*p10)(mytimes*4, myfit10 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot10 <- plotCombined(mypred10, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot10)
# myplot10
##   index     value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1    72 -2356.519      TRUE         99         -5.984038     -3.407580
## 2    23 -2356.519      TRUE        100         -5.984055     -3.407593
## 3    90 -2356.519      TRUE         98         -5.984039     -3.407588
## 4   144 -2356.519     FALSE        100         -5.984040     -3.407582
## 5    63 -2356.519      TRUE         78         -5.983866     -3.407695
## 6   135 -2356.519     FALSE        100         -5.984081     -3.407616
##   logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1         -5.782424         -10.57408     -6.696947  -7.013798
## 2         -5.782408         -10.57408     -6.696946  -7.013783
## 3         -5.782429         -10.57408     -6.696951  -7.013733
## 4         -5.782424         -10.57408     -6.696946  -7.013786
## 5         -5.782736         -10.57412     -6.697111  -7.013753
## 6         -5.782389         -10.57409     -6.696947  -7.013793
##   log_srel_apap_Chan1997_14_0 log_s0_apap_Chan1997_14_0
## 1                    3.289950                 -7.052483
## 2                    3.289934                 -7.052496
## 3                    3.289921                 -7.052506
## 4                    3.289937                 -7.052506
## 5                    3.289810                 -7.052222
## 6                    3.289991                 -7.052427
##   log_srel_apap_glu_Chan1997_14_0 log_s0_apap_glu_Chan1997_14_0
## 1                        4.404972                     -7.517436
## 2                        4.405154                     -7.517363
## 3                        4.405028                     -7.517374
## 4                        4.405060                     -7.517341
## 5                        4.403288                     -7.516821
## 6                        4.405348                     -7.517387
##   log_srel_apap_sul_Chan1997_14_0 log_s0_apap_sul_Chan1997_14_0
## 1                        4.959255                     -8.122914
## 2                        4.958957                     -8.123157
## 3                        4.959125                     -8.123015
## 4                        4.959096                     -8.123068
## 5                        4.960016                     -8.122510
## 6                        4.958673                     -8.123290
##   log_srel_apap_cys_Chan1997_14_0 log_s0_apap_cys_Chan1997_14_0
## 1                       -6.968630                     -9.728692
## 2                       -3.817880                     -9.728688
## 3                       -5.081099                     -9.728682
## 4                       -3.199117                     -9.728689
## 5                       -2.727214                     -9.728656
## 6                       -2.482580                     -9.728732
##   log_srel_apap_Chiew2010_56_0 log_s0_apap_Chiew2010_56_0
## 1                    0.1017979                  -7.141670
## 2                    0.1020375                  -7.141873
## 3                    0.1016049                  -7.141541
## 4                    0.1018871                  -7.141784
## 5                    0.1016913                  -7.141588
## 6                    0.1019149                  -7.141807
##   log_srel_apap_glu_Chiew2010_56_0 log_s0_apap_glu_Chiew2010_56_0
## 1                         5.497352                      -6.479678
## 2                         5.497414                      -6.479689
## 3                         5.497194                      -6.479812
## 4                         5.497293                      -6.479720
## 5                         5.496531                      -6.479613
## 6                         5.498219                      -6.479841
##   log_srel_apap_sul_Chiew2010_56_0 log_s0_apap_sul_Chiew2010_56_0
## 1                         2.464006                      -6.105482
## 2                         2.463693                      -6.105647
## 3                         2.463902                      -6.105591
## 4                         2.463910                      -6.105581
## 5                         2.468499                      -6.106767
## 6                         2.463066                      -6.105751
##   log_srel_apap_Critchley2005_14_0 log_s0_apap_Critchley2005_14_0
## 1                         2.350997                      -6.825250
## 2                         2.350973                      -6.825248
## 3                         2.350992                      -6.825255
## 4                         2.350957                      -6.825260
## 5                         2.350889                      -6.824922
## 6                         2.350855                      -6.825138
##   log_srel_apap_glu_Critchley2005_14_0 log_s0_apap_glu_Critchley2005_14_0
## 1                             3.602125                          -6.967212
## 2                             3.602291                          -6.967185
## 3                             3.602170                          -6.967135
## 4                             3.602175                          -6.967151
## 5                             3.601192                          -6.966273
## 6                             3.602554                          -6.967230
##   log_srel_apap_sul_Critchley2005_14_0 log_s0_apap_sul_Critchley2005_14_0
## 1                             4.017564                          -7.716869
## 2                             4.017535                          -7.716892
## 3                             4.017606                          -7.716842
## 4                             4.017561                          -7.716888
## 5                             4.017116                          -7.717051
## 6                             4.017688                          -7.716863
##   log_srel_apap_cys_Critchley2005_14_0 log_s0_apap_cys_Critchley2005_14_0
## 1                            -8.051930                          -10.10211
## 2                            -4.902004                          -10.10212
## 3                            -6.164908                          -10.10211
## 4                            -4.283094                          -10.10212
## 5                            -3.755340                          -10.10221
## 6                            -3.566113                          -10.10215
##   log_srel_apap_Rawlins1977_0_1 log_s0_apap_Rawlins1977_0_1
## 1                      2.541371                   -8.725966
## 2                      2.541383                   -8.725981
## 3                      2.541391                   -8.725948
## 4                      2.541383                   -8.725955
## 5                      2.541252                   -8.725793
## 6                      2.541297                   -8.725894
##   log_srel_apap_Rawlins1977_1_0 log_s0_apap_Rawlins1977_1_0
## 1                      2.772759                   -8.547667
## 2                      2.772043                   -8.546641
## 3                      2.770963                   -8.545399
## 4                      2.772593                   -8.547197
## 5                      2.772166                   -8.547219
## 6                      2.772761                   -8.547107
##   log_srel_apap_Rawlins1977_2_0 log_s0_apap_Rawlins1977_2_0
## 1                      2.485538                   -8.345503
## 2                      2.485019                   -8.343667
## 3                      2.485324                   -8.344722
## 4                      2.485179                   -8.344683
## 5                      2.485390                   -8.344672
## 6                      2.485404                   -8.344505
##   log_srel_apap_Rawlins1977_05_0 log_s0_apap_Rawlins1977_05_0
## 1                       3.075085                    -8.175010
## 2                       3.075391                    -8.175087
## 3                       3.075176                    -8.175017
## 4                       3.074992                    -8.174969
## 5                       3.074871                    -8.175119
## 6                       3.074939                    -8.175055
##                    logAPAPCYS_HLM_CL:  2.5570243395892e-05
##                        logAPAPCYS_Km:  0.00123467580981531
##                    logAPAPGLU_HLM_CL:  0.00251863572465911
##                        logAPAPGLU_Km:  0.033121273391679
##                    logAPAPSUL_HLM_CL:  0.00308123765707667
##                           logKa_apap:  0.000899386426280096
##            log_s0_apap_Chan1997_14_0:  0.000865258189660525
##           log_s0_apap_Chiew2010_56_0:  0.000791429248666158
##       log_s0_apap_Critchley2005_14_0:  0.00108600472415769
##        log_s0_apap_cys_Chan1997_14_0:  5.9550106934986e-05
##   log_s0_apap_cys_Critchley2005_14_0:  4.09928708157083e-05
##        log_s0_apap_glu_Chan1997_14_0:  0.000543524495738245
##       log_s0_apap_glu_Chiew2010_56_0:  0.00153430475174298
##   log_s0_apap_glu_Critchley2005_14_0:  0.000942275918358174
##          log_s0_apap_Rawlins1977_0_1:  0.000162315884677548
##         log_s0_apap_Rawlins1977_05_0:  0.000281603560301461
##          log_s0_apap_Rawlins1977_1_0:  0.000193997116676406
##          log_s0_apap_Rawlins1977_2_0:  0.000237462036719598
##        log_s0_apap_sul_Chan1997_14_0:  0.000296662808440626
##       log_s0_apap_sul_Chiew2010_56_0:  0.00223060565323054
##   log_s0_apap_sul_Critchley2005_14_0:  0.000445252327814248
##          log_srel_apap_Chan1997_14_0:  26.8415245012767
##         log_srel_apap_Chiew2010_56_0:  1.10715971018415
##     log_srel_apap_Critchley2005_14_0:  10.4960243104907
##      log_srel_apap_cys_Chan1997_14_0:  0.00094094138090811
## log_srel_apap_cys_Critchley2005_14_0:  0.000318486578636513
##      log_srel_apap_glu_Chan1997_14_0:  81.8568173051488
##     log_srel_apap_glu_Chiew2010_56_0:  244.044756600744
## log_srel_apap_glu_Critchley2005_14_0:  36.6760946617877
##        log_srel_apap_Rawlins1977_0_1:  12.6970698875144
##       log_srel_apap_Rawlins1977_05_0:  21.6517119661049
##        log_srel_apap_Rawlins1977_1_0:  16.0027176603338
##        log_srel_apap_Rawlins1977_2_0:  12.0075810603921
##      log_srel_apap_sul_Chan1997_14_0:  142.487667782381
##     log_srel_apap_sul_Chiew2010_56_0:  11.7517960003797
## log_srel_apap_sul_Critchley2005_14_0:  55.5656002676397

Data revisited

Some sigmas of the data are NA, try to recover an estimate for the sigmas from the other sigmas.

Take out some outliers for the fitting. Those are: study, name, time, reason 1. Chan1997, apap_cys, 5400, don’t know what went wrong with this one, but it just doesn’t fit reasonably in the time course 3. Chan1997, apap, 18000, The sigma is a few orders of magnitude lower. Mirjam said it might be that in this point less people were measured (eg 2) and they had nearly the same value. Then of course, sigma would be very smal

# data <- data %>% 
#   filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>% 
#   filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>% 
#   {.}

 data %>% select(-n) %>% filter(study %>% str_detect("Chiew"))%>% as.datalist() %>% plotData()

data_with_errors <- data %>% 
  filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>% 
  filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>% 
  fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name"), blather = T)

data_with_errors %>% select(study, name, D_apap, s0, srel) %>% unique()
##             study     name D_apap         s0       srel
## 1        Chan1997     apap    1.4 -14.395145  -4.516431
## 13       Chan1997 apap_sul    1.4 -15.377176  -2.609451
## 26       Chan1997 apap_glu    1.4 -14.970003  -1.709818
## 39       Chan1997 apap_cys    1.4 -19.186945 -13.879472
## 51      Chiew2010     apap    5.6 -24.503629  -1.232431
## 65      Chiew2010 apap_glu    5.6  -9.710964  -2.303673
## 79      Chiew2010 apap_sul    5.6 -10.576758  -4.220755
## 93  Critchley2005     apap    1.4 -11.298969 -16.644666
## 108 Critchley2005 apap_sul    1.4 -16.328137  -2.802784
## 123 Critchley2005 apap_glu    1.4 -15.826025  -2.457868
## 138 Critchley2005 apap_cys    1.4 -18.851664 -18.687321
## 153   Rawlins1977     apap      0 -13.826728  -2.924789
## 163   Rawlins1977     apap      1 -28.577775  -1.968079
## 171   Rawlins1977     apap      2 -27.097492  -1.002046
## 179   Rawlins1977     apap    0.5 -27.374147  -2.115781

What’s the problem with the Chiew-dataset? The error model doesn’t work out, somehow

data_with_errors %>% 
  filter(study %>% str_detect("Chiew"), name %in% "apap") %>% 
ggplot(aes(x=value)) +
            geom_point(aes(y=sigmaLS^2*(n), color = log(time))) +
            geom_line(aes(y=sigma^2*n)) +
            geom_ribbon(aes(ymin=cbLower95, ymax=cbUpper95), alpha=.3) +
            geom_ribbon(aes(ymin=cbLower68, ymax=cbUpper68), alpha=.3) +
            ylab("variance") +
            facet_wrap(~condidnt, scales = "free") +
            scale_y_log10() +
            theme_dMod() + 
            scale_color_continuous( low = "#98f5ff", high = "#4c4cdb")
myplot <- data_with_errors %>% 
  rename(sigma_fitted = sigma) %>% 
  left_join(data) %>% 
  gather("which_sig", "sigma", sigma, sigma_fitted) %>% 
  mutate(condition = paste0(study, D_apap)) %>% 
  ggplot(aes(x = time,y = log10(sigma))) +
  geom_line(aes(color = condition, linetype = which_sig)) +
  # geom_point(aes(color = condition, shape = which_sig))+
  facet_wrap("name", scales = "free") 
## Joining, by = c("study", "n", "time", "name", "value", "D_apap", "Ave_apap")
myplot %>% plotly::ggplotly()
# myplot

Fit error model first, then fit the model with free scaling factors

load("methacetin.rda")

x <- Xs(myodemodel) # make prediction function
loadDLL(x)

# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]

free_parameters_e1 <- c("APAPGLU_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPSUL_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPCYS_HLM_CL",  # Vmax value
                     "APAPCYS_Km",  # Km value

                     "Ka_apap"#, #"F_apap_sul"    ,
                     # "Kpre_apap", "Kpki_apap", "Kpli_apap",
                     # "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
                     # "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
                     # "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
                     # "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
                     # "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
                     )

fixed_parameters_e1 <- pars[!(names(pars)%in%c(free_parameters_e1,names(f)[1]))] %>% names

mydatalist <- data %>% 
  filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>% 
  filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>% 
  fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name")) %>% 
  select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")


observables_e1 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
                 apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
                 apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
                 apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_e1 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)

# free_parameters_e1 <- c(free_parameters_e1, scale_parameters_e1)

# error_model_e1 <- c(apap = "srel_apap*apap^2 +s0_apap", 
#                  apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu", 
#                  apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul", 
#                  apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_e1 <- setdiff(getSymbols(error_model_e1), names(error_model_e1)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
  cond <- unlist(conditions[i,])[2:3]
  
  trafo <- as.character(pars) %>% set_names(names(pars))
  trafo[names(cond)] <- cond
  trafo[free_parameters_e1] <- paste0("exp(log", free_parameters_e1, ")")

  scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_e1, x = scale_parameters_e1, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_e1)
  scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  
  # errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_e1, x = error_parameters_e1, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_e1)
  # errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  
  trafo <- c(trafo, scales)#, errors)
    
  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p_e1 <- NULL
for(i in 1:length(p_list)) { p_e1 <<- p_e1 + p_list[[i]]}

g_e1 <- Y(observables_e1, x)#, parameters = c(free_parameters_e1, scale_parameters_e1))  

# err_e1 <- Y(error_model_e1, g_e1)

obj_e1 <- normL2(mydatalist, (g_e1*x*p_e1))#, errmodel = err_e1)


pouter_e1 <- rep(0, length(getParameters(obj_e1))) %>% set_names(getParameters(obj_e1))
# pouter_e1[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()

#job_e1 <- runbg({myfit <- mstrust(objfun = obj_e1, center = pouter_e1, sd = 3, studyname = "methacetin", cores = 24, fits = 200); myfit}, machine = "ruprecht1", filename = "job_e1", input = global_env_without(c("mypred", "job", "myfit")))

# save(job_e1, file = "job_e1.rda")

# job_e1$check()
# load("job_e1_1_result.RData")
# myfit_e1 <- .runbgOutput
# load("job_e1.rda")
# myfit_e1 <- job_e1$get()
# save(myfit_e1, file = "myfit_e1.rda")
# job_e1$purge()
load("myfit_e1.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit_e1 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_e1 %>% as.parframe())#+scale_y_log10()

mypred_e1 <- (g_e1*x*p_e1)(mytimes*4, myfit_e1 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot_e1 <- plotCombined(mypred_e1, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_e1)
# myplot_e1
##   index    value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1   135 1151.996      TRUE         69        -18.255493     -2.247950
## 2    70 1151.996     FALSE        100        -25.864555     -2.246878
## 3   102 1151.996      TRUE         87        -20.902541     -2.247412
## 4   171 1151.996      TRUE         61        -15.814734     -2.245022
## 5   191 1152.002      TRUE         60         -7.374898     -2.255883
## 6    81 1152.021     FALSE        100         -5.978911     -2.280757
##   logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1         -4.439263         -9.981200     -7.521667  -7.664488
## 2         -4.439394         -9.980696     -7.521636  -7.664391
## 3         -4.439279         -9.981359     -7.521708  -7.664394
## 4         -4.439563         -9.979954     -7.521437  -7.664217
## 5         -4.482247         -9.973992     -7.522330  -7.663926
## 6         -4.629608         -9.955660     -7.524306  -7.662562
##   log_scale_apap_Chan1997_14_0 log_scale_apap_glu_Chan1997_14_0
## 1                    0.6599173                       12.8364130
## 2                    0.6598512                       20.4463668
## 3                    0.6599152                       15.4839438
## 4                    0.6597140                       10.3980699
## 5                    0.6598690                        1.9483154
## 6                    0.6599017                        0.5291009
##   log_scale_apap_sul_Chan1997_14_0 log_scale_apap_cys_Chan1997_14_0
## 1                       -1.0771287                       -0.4991120
## 2                       -1.0770584                       -0.4996198
## 3                       -1.0770897                       -0.4989502
## 4                       -1.0770215                       -0.5003547
## 5                       -1.0342097                       -0.5064348
## 6                       -0.8868909                       -0.5250863
##   log_scale_apap_Chiew2010_56_0 log_scale_apap_glu_Chiew2010_56_0
## 1                     0.4747969                         14.360289
## 2                     0.4746995                         21.970031
## 3                     0.4747278                         17.007724
## 4                     0.4745103                         11.921363
## 5                     0.4764708                          3.475552
## 6                     0.4813080                          2.066233
##   log_scale_apap_sul_Chiew2010_56_0 log_scale_apap_Critchley2005_14_0
## 1                        -0.8302279                         0.4550537
## 2                        -0.8301517                         0.4549816
## 3                        -0.8301895                         0.4550346
## 4                        -0.8301000                         0.4548393
## 5                        -0.7857834                         0.4549590
## 6                        -0.6342932                         0.4548708
##   log_scale_apap_glu_Critchley2005_14_0
## 1                            13.2497621
## 2                            20.8597371
## 3                            15.8973054
## 4                            10.8114742
## 5                             2.3620052
## 6                             0.9437402
##   log_scale_apap_sul_Critchley2005_14_0
## 1                             -1.373527
## 2                             -1.373459
## 3                             -1.373493
## 4                             -1.373425
## 5                             -1.330606
## 6                             -1.183280
##   log_scale_apap_cys_Critchley2005_14_0 log_scale_apap_Rawlins1977_0_1
## 1                            -0.7385315                      0.1580270
## 2                            -0.7390379                      0.1579934
## 3                            -0.7383715                      0.1580312
## 4                            -0.7397729                      0.1579200
## 5                            -0.7458288                      0.1579958
## 6                            -0.7644113                      0.1579937
##   log_scale_apap_Rawlins1977_1_0 log_scale_apap_Rawlins1977_2_0
## 1                      0.3682622                      0.1562957
## 2                      0.3682065                      0.1562248
## 3                      0.3682818                      0.1562731
## 4                      0.3680777                      0.1560860
## 5                      0.3680259                      0.1564571
## 6                      0.3675088                      0.1570966
##   log_scale_apap_Rawlins1977_05_0
## 1                       0.1238156
## 2                       0.1237872
## 3                       0.1238273
## 4                       0.1236985
## 5                       0.1235803
## 6                       0.1229841
##                     logAPAPCYS_HLM_CL:  4.6261539108133e-05
##                         logAPAPCYS_Km:  0.000541229613400164
##                     logAPAPGLU_HLM_CL:  1.17961442872092e-08
##                         logAPAPGLU_Km:  0.105615516447963
##                     logAPAPSUL_HLM_CL:  0.0118046378599577
##                            logKa_apap:  0.000469196966461811
##          log_scale_apap_Chan1997_14_0:  1.93463224874753
##         log_scale_apap_Chiew2010_56_0:  1.60768769704826
##     log_scale_apap_Critchley2005_14_0:  1.57625801184958
##      log_scale_apap_cys_Chan1997_14_0:  0.607069478303199
## log_scale_apap_cys_Critchley2005_14_0:  0.47781508039879
##      log_scale_apap_glu_Chan1997_14_0:  375649.964718959
##     log_scale_apap_glu_Chiew2010_56_0:  1724226.50632285
## log_scale_apap_glu_Critchley2005_14_0:  567934.927599969
##        log_scale_apap_Rawlins1977_0_1:  1.1711977984144
##       log_scale_apap_Rawlins1977_05_0:  1.13180719602812
##        log_scale_apap_Rawlins1977_1_0:  1.44522093867267
##        log_scale_apap_Rawlins1977_2_0:  1.16917192088704
##      log_scale_apap_sul_Chan1997_14_0:  0.340572019805132
##     log_scale_apap_sul_Chiew2010_56_0:  0.435949924368047
## log_scale_apap_sul_Critchley2005_14_0:  0.253212324969791

Fit error model to data, then fit model without scaling

load("methacetin.rda")

x <- Xs(myodemodel) # make prediction function
loadDLL(x)

# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]

free_parameters_e2 <- c("APAPGLU_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPSUL_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPCYS_HLM_CL",  # Vmax value
                     "APAPCYS_Km",  # Km value

                     "Ka_apap"#, #"F_apap_sul"    ,
                     # "Kpre_apap", "Kpki_apap", "Kpli_apap",
                     # "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
                     # "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
                     # "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
                     # "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
                     # "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
                     )

fixed_parameters_e2 <- pars[!(names(pars)%in%c(free_parameters_e2,names(f)[1]))] %>% names

mydatalist <- data %>% 
  filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>% 
  filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>% 
  fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name")) %>% 
  select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")


observables_e2 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
                 apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
                 apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
                 apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_e2 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)

# free_parameters_e2 <- c(free_parameters_e2, scale_parameters_e2)

# error_model_e2 <- c(apap = "srel_apap*apap^2 +s0_apap", 
#                  apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu", 
#                  apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul", 
#                  apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_e2 <- setdiff(getSymbols(error_model_e2), names(error_model_e2)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
  cond <- unlist(conditions[i,])[2:3]
  
  trafo <- as.character(pars) %>% set_names(names(pars))
  trafo[names(cond)] <- cond
  trafo[free_parameters_e2] <- paste0("exp(log", free_parameters_e2, ")")

  scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_e2, x = scale_parameters_e2, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_e2)
  # scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  scales[1:length(scales)] <- "1"
  
  # errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_e2, x = error_parameters_e2, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_e2)
  # errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  
  trafo <- c(trafo, scales)#, errors)
    
  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p_e2 <- NULL
for(i in 1:length(p_list)) { p_e2 <<- p_e2 + p_list[[i]]}

g_e2 <- Y(observables_e2, x)#, parameters = c(free_parameters_e2, scale_parameters_e2))  

# err_e2 <- Y(error_model_e2, g_e2)

obj_e2 <- normL2(mydatalist, (g_e2*x*p_e2))#, errmodel = err_e2)


pouter_e2 <- rep(0, length(getParameters(obj_e2))) %>% set_names(getParameters(obj_e2))
# pouter_e2[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()

# job_e2 <- runbg({myfit <- mstrust(objfun = obj_e2, center = pouter_e2, sd = 5, studyname = "methacetin", cores = 20, fits = 400); myfit}, machine = "ruprecht2", filename = "job_e2", input = global_env_without(c("mypred", "job", "myfit", "myplot")))

# save(job_e2, file = "job_e2.rda")

# job_e2$check()
# load("job_e2_1_result.RData")
# myfit_e2 <- .runbgOutput
# load("job_e1.rda")
# myfit_e1 <- job_e1$get()
# save(myfit_e2, file = "myfit_e2.rda")
# job_e1$purge()
load("myfit_e2.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit_e2 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_e2 %>% as.parframe())#+scale_y_log10()

mypred_e2 <- (g_e2*x*p_e2)(mytimes*4, myfit_e2 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot_e2 <- plotCombined(mypred_e2, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_e2)
# myplot_e2
##   index    value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1   379 2402.825      TRUE         56         -6.195743     -3.728993
## 2   125 2402.825      TRUE         46         -6.195869     -3.729162
## 3   147 2402.826      TRUE         68         -6.195203     -3.728064
## 4   224 2402.826      TRUE         54         -6.195879     -3.729148
## 5   374 2402.826      TRUE         39         -6.195582     -3.728692
## 6    28 2402.826      TRUE         38         -6.195853     -3.729130
##   logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1         -5.922382         -10.65531     -7.895205  -7.156523
## 2         -5.922387         -10.65532     -7.895313  -7.156831
## 3         -5.922448         -10.65529     -7.894827  -7.156478
## 4         -5.922413         -10.65534     -7.895534  -7.156525
## 5         -5.922399         -10.65530     -7.895033  -7.156939
## 6         -5.922391         -10.65534     -7.895523  -7.156850
## logAPAPCYS_HLM_CL:  2.3575301587201e-05
##     logAPAPCYS_Km:  0.000372525594735151
## logAPAPGLU_HLM_CL:  0.00203808917141075
##     logAPAPGLU_Km:  0.0240169974041635
## logAPAPSUL_HLM_CL:  0.00267881249032963
##        logKa_apap:  0.000779761067836627

Tackle too slow dynamics of apap_*

The dynamics of e.g. apap_glu are too slow. Test two hypotheses:

  • Increasing CLrenal_apap* fixes it
  • multiplying a rate constant to Ave_apap_ -> A_lu_apap_ fixes it.

Compile a model with less sensitivities for faster integration

This model also comprises some rate constants to make some reactions go faster

# Multiply rates from "ve" to "lu" with additional free parameters
reactions_to_lung <- reactions
reactions_to_lung[reactions_to_lung$from %>% str_detect("(1*Ave_apap)"),3] %<>% paste0(" * (k_to_lung_", c("sul", "", "cys", "glu"),")") 

reactions_to_lung
el_to_lung <- eqnlist()

for(i in 1:nrow(reactions_to_lung)) el_to_lung <- addReaction(el_to_lung, reactions_to_lung$from[i], reactions_to_lung$to[i], reactions_to_lung$rate[i], reactions_to_lung$description[i])


# Convert to "eqnvec", which is basically a named vector of the ODEs and the names denote the states
f_to_lung <- el_to_lung %>% as.eqnvec()

parameters <- all_pars[getParameters(x)] %>% names()
free_parameters_to_lung <- parameters[parameters %>% sapply(. %>% str_detect(c("CL", "Km", "Ka_apap$")) %>% any) ]

# myodemodel_to_lung <- odemodel(f_to_lung, modelname = "methacetin_to_lung", fixed = setdiff(parameters, c(free_parameters_to_lung, "Ave_apap"))) # Include "Ave_apap" because of a bug: You need the sensitivities to at least one initial value
save(myodemodel_to_lung, file = "methacetin_to_lung.rda")

Fit the model with CLrenal*-parameters The hypothesis was that the slowly decaying dynamics of _glu, _cys and _sul could be fixed by removing them faster from the system due to clearance by the kidney. This didn’t work.

load("methacetin_to_lung.rda")

x_to_lung <- Xs(myodemodel_to_lung) # make prediction function
loadDLL(x_to_lung)

# get the only the parameters needed for x
pars_to_lung <- rep(1,4) %>% set_names(paste0("k_to_lung_", c("sul", "", "cys", "glu")))

pars <- c(all_pars, pars_to_lung)[getParameters(x_to_lung)]

free_parameters_CLrenal <- c("APAPGLU_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPSUL_HLM_CL",  # Vmax value
                     "APAPGLU_Km",  # Km value
                     "APAPCYS_HLM_CL",  # Vmax value
                     "APAPCYS_Km",  # Km value

                     "Ka_apap",
                     "CLrenal_apap", "CLrenal_apap_cys", "CLrenal_apap_glu", "CLrenal_apap_sul"
                     )

fixed_parameters_CLrenal <- pars[!(names(pars)%in%c(free_parameters_CLrenal,names(f)[1]))] %>% names

mydatalist <- data %>% 
  filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>% 
  filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>% 
  fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name")) %>% 
  select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")


observables_CLrenal <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
                 apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
                 apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
                 apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_CLrenal <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)

# free_parameters_CLrenal <- c(free_parameters_CLrenal, scale_parameters_CLrenal)

# error_model_CLrenal <- c(apap = "srel_apap*apap^2 +s0_apap", 
#                  apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu", 
#                  apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul", 
#                  apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_CLrenal <- setdiff(getSymbols(error_model_CLrenal), names(error_model_CLrenal)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
  cond <- unlist(conditions[i,])[2:3]
  
  trafo <- as.character(pars) %>% set_names(names(pars))
  trafo[names(cond)] <- cond
  trafo[free_parameters_CLrenal] <- paste0("exp(log", free_parameters_CLrenal, ")")

  scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_CLrenal, x = scale_parameters_CLrenal, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_CLrenal)
  # scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  scales[1:length(scales)] <- "1"
  
  # errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_CLrenal, x = error_parameters_CLrenal, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_CLrenal)
  # errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  
  trafo <- c(trafo, scales)#, errors)
    
  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p_CLrenal <- NULL
for(i in 1:length(p_list)) { p_CLrenal <<- p_CLrenal + p_list[[i]]}

g_CLrenal <- Y(observables_CLrenal, x_to_lung)#, parameters = c(free_parameters_CLrenal, scale_parameters_CLrenal))  

# err_CLrenal <- Y(error_model_CLrenal, g_CLrenal)

obj_CLrenal <- normL2(mydatalist, (g_CLrenal*x_to_lung*p_CLrenal))#, errmodel = err_CLrenal)


pouter_CLrenal <- rep(0, length(getParameters(obj_CLrenal))) %>% set_names(getParameters(obj_CLrenal))
pouter_CLrenal[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
pouter_CLrenal[names(pouter_CLrenal) %>% str_detect("CLr")] <- -4
pouter_CLrenal[1:2] <- c(1.5,5)
# obj_CLrenal(pouter_CLrenal)
# job_CLrenal <- runbg({myfit <- mstrust(objfun = obj_CLrenal, center = pouter_CLrenal, sd = 3, studyname = "methacetin", cores = 16, fits = 100, iterlim = 150); myfit}, machine = "ruprecht2", filename = "job_CLrenal", input = global_env_without(c("mypred", "job", "myfit", "myplot")))

# save(job_CLrenal, file = "job_CLrenal.rda")

# job_CLrenal$check()
# myfit_CLrenal <- job_CLrenal$get()$ruprecht2
# # save(myfit_CLrenal, file = "myfit_CLrenal.rda")
# # job_CLrenal$purge()
load("myfit_CLrenal.rda")
loadDLL(x_to_lung)
## The following local files were dynamically loaded: methacetin_to_lung.so, methacetin_to_lung_s.so
myfit_CLrenal %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_CLrenal %>% as.parframe())+scale_y_log10()

mypouter <- myfit_CLrenal %>% as.parframe() %>% as.parvec
# mypouter[1:2] <- mypouter[1:2]+1
mypred_CLrenal <- (g_CLrenal*x_to_lung*p_CLrenal)(mytimes*4, mypouter, deriv = F)

# myplot_CLrenal <- plotCombined(mypred_CLrenal, mydatalist, name %>% sapply(. %>% str_detect(c(paste0("^",names(observables)), "Agu", "Ali", "Ave")) %>% any))
myplot_CLrenal <- plotCombined(mypred_CLrenal, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_CLrenal)
# myplot_CLrenal
##   index    value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1    22 1679.527      TRUE         55         -6.532773     -4.069886
## 2     7 1679.527      TRUE         46         -6.533129     -4.070574
## 3    60 1679.527      TRUE         59         -6.532886     -4.069996
## 4    31 1679.527      TRUE         60         -6.532866     -4.069961
## 5    70 1679.527      TRUE         55         -6.532908     -4.070019
## 6    71 1679.527      TRUE         66         -6.532916     -4.070109
##   logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logCLrenal_apap_sul
## 1         -5.629237         -10.64046     -6.378762           -5.099622
## 2         -5.629241         -10.64043     -6.374378           -5.099647
## 3         -5.629227         -10.64037     -6.382149           -5.099616
## 4         -5.629244         -10.64037     -6.383058           -5.099634
## 5         -5.629232         -10.64037     -6.383474           -5.099626
## 6         -5.629243         -10.64041     -6.380751           -5.099644
##   logKa_apap logCLrenal_apap logCLrenal_apap_cys logCLrenal_apap_glu
## 1  -7.220531       -32.01661           -5.271552           -5.448181
## 2  -7.220536       -28.60979           -5.273667           -5.448045
## 3  -7.220632       -26.47175           -5.269753           -5.448253
## 4  -7.220556       -33.60964           -5.269310           -5.448251
## 5  -7.220625       -19.92413           -5.269103           -5.448267
## 6  -7.220526       -27.18079           -5.270518           -5.448182
##   logAPAPCYS_HLM_CL:  2.39279181355532e-05
##       logAPAPCYS_Km:  0.0016972225659718
##   logAPAPGLU_HLM_CL:  0.00145496514549494
##       logAPAPGLU_Km:  0.0170793409249032
##   logAPAPSUL_HLM_CL:  0.00359131615388076
##     logCLrenal_apap:  1.24556095423351e-14
## logCLrenal_apap_cys:  0.00513563209480641
## logCLrenal_apap_glu:  0.0043041255170805
## logCLrenal_apap_sul:  0.00609905356785048
##          logKa_apap:  0.000731413761365896

Allow rate from venous to lung to be faster - 2 local optima in the plot

Chiew2010 apap_sul looks good for the first time, even without scaling

load("methacetin_to_lung.rda")

x_to_lung <- Xs(myodemodel_to_lung) # make prediction function
loadDLL(x_to_lung)

# get the only the parameters needed for x
pars_to_lung <- rep(0,4) %>% set_names(paste0("k_to_lung_", c("sul", "", "cys", "glu")))

pars <- c(all_pars, pars_to_lung)[getParameters(x)]

free_parameters_to_lung <- c("APAPGLU_HLM_CL", "APAPGLU_Km", "APAPSUL_HLM_CL", "APAPSUL_Km", "APAPCYS_HLM_CL", "APAPCYS_Km", 
                             "Ka_apap", 
                             paste0("k_to_lung_", c("sul", "", "cys", "glu"))
                             # "CLrenal_apap", "CLrenal_metc13",  "CLrenal_apap_sul", "CLrenal_apap_cys", "CLrenal_apap_glu", "CLrenal_co2c13"
                             )


mydatalist <- data %>% 
  filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>% 
  filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>% 
  fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name")) %>% 
  select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")


observables_to_lung <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
                 apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
                 apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
                 apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_to_lung <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)

# free_parameters_to_lung <- c(free_parameters_to_lung, scale_parameters_to_lung)

# error_model_to_lung <- c(apap = "srel_apap*apap^2 +s0_apap", 
#                  apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu", 
#                  apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul", 
#                  apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_to_lung <- setdiff(getSymbols(error_model_to_lung), names(error_model_to_lung)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
  cond <- unlist(conditions[i,])[2:3]
  
  trafo <- as.character(pars) %>% set_names(names(pars))
  trafo[names(cond)] <- cond
  trafo[free_parameters_to_lung] <- paste0("exp(log", free_parameters_to_lung, ")")

  scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_to_lung, x = scale_parameters_to_lung, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_to_lung)
  # scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  scales[1:length(scales)] <- "1"
  
  # errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_to_lung, x = error_parameters_to_lung, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_to_lung)
  # errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  
  trafo <- c(trafo, scales)#, errors)
    
  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p_to_lung <- NULL
for(i in 1:length(p_list)) { p_to_lung <<- p_to_lung + p_list[[i]]}

g_to_lung <- Y(observables_to_lung, x_to_lung)#, parameters = c(free_parameters_to_lung, scale_parameters_to_lung))  

# err_to_lung <- Y(error_model_to_lung, g_to_lung)

obj_to_lung <- normL2(mydatalist, (g_to_lung*x_to_lung*p_to_lung))#, errmodel = err_to_lung)


pouter_to_lung <- rep(0, length(getParameters(obj_to_lung))) %>% set_names(getParameters(obj_to_lung))
pouter_to_lung[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
obj_to_lung(pouter_to_lung)
# job_to_lung <- runbg({myfit <- mstrust(objfun = obj_to_lung, center = pouter_to_lung, sd = 2, studyname = "methacetin", cores = 24, fits = 100, iterlim = 150); myfit}, machine = "ruprecht1", filename = "job_to_lung", input = global_env_without(c("mypred", "job", "myfit", "myplot")))

# save(job_to_lung, file = "job_to_lung.rda")

# job_to_lung$check()
# myfit_to_lung <- job_to_lung$get()$ruprecht1
# save(myfit_to_lung, file = "myfit_to_lung.rda")
# job_to_lung$purge()
load("myfit_to_lung.rda")
loadDLL(x_to_lung)
## The following local files were dynamically loaded: methacetin_to_lung.so, methacetin_to_lung_s.so
myfit_to_lung %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_to_lung %>% as.parframe())+scale_y_log10()

mypred_to_lung <- (g_to_lung*x_to_lung*p_to_lung)(mytimes*4, myfit_to_lung %>% as.parframe() %>% as.parvec(1) %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)

mypred_to_lung2 <- (g_to_lung*x_to_lung*p_to_lung)(mytimes*4, myfit_to_lung %>% as.parframe() %>% as.parvec(55) %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F) %>% set_names(paste(names(.),"_2"))
# myplot_to_lung <- plotCombined(mypred_to_lung, mydatalist, name %>% sapply(. %>% str_detect(c(paste0("^",names(observables)), "Agu", "Ali", "Ave")) %>% any))

myplot_to_lung <- plotCombined(c(mypred_to_lung,mypred_to_lung2), mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_to_lung)
# myplot_to_lung
##   index    value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1    70 1557.189      TRUE        142          1.882061      5.800306
## 2    34 1557.199      TRUE        148          1.296540      5.215847
## 3    41 1557.213      TRUE        142          1.652053      5.570782
## 4     4 1557.214      TRUE        150          1.600278      5.520080
## 5    54 1557.249     FALSE        150          1.315475      5.237448
## 6    95 1557.256     FALSE        150          1.215461      5.132848
##   logAPAPSUL_HLM_CL logAPAPSUL_Km logAPAPCYS_HLM_CL logAPAPCYS_Km
## 1         1.3139749      5.074520         -7.348043     -6.117347
## 2         2.0624390      5.822052         -7.348164     -6.117514
## 3         1.0046156      4.764574         -7.348081     -6.117433
## 4         1.0356908      4.794728         -7.348121     -6.117488
## 5         0.7347650      4.491662         -7.348263     -6.117822
## 6         0.7286464      4.489584         -7.348001     -6.117180
##   logk_to_lung_sul logKa_apap logk_to_lung_ logk_to_lung_cys
## 1       -0.1459383  -7.460770    -0.2039797         3.060832
## 2       -0.1448628  -7.460689    -0.2039490         3.060762
## 3       -0.1454166  -7.460765    -0.2039929         3.060800
## 4       -0.1444364  -7.460705    -0.2039785         3.060779
## 5       -0.1423281  -7.460575    -0.2039709         3.060692
## 6       -0.1464887  -7.460877    -0.2040240         3.060828
##   logk_to_lung_glu
## 1       -0.8190345
## 2       -0.8201550
## 3       -0.8195732
## 4       -0.8207004
## 5       -0.8229671
## 6       -0.8182854
## logAPAPCYS_HLM_CL:  0.00064385123805241
##     logAPAPCYS_Km:  0.0022042952730362
## logAPAPGLU_HLM_CL:  6.56702557689903
##     logAPAPGLU_Km:  330.40080787627
## logAPAPSUL_HLM_CL:  3.72093453867929
##     logAPAPSUL_Km:  159.895379325816
##        logKa_apap:  0.000575213236929382
##     logk_to_lung_:  0.815478948433962
##  logk_to_lung_cys:  21.3453185074998
##  logk_to_lung_glu:  0.44085707897448
##  logk_to_lung_sul:  0.864211035946343

increase rate to lung and allow scaling - didn’t work out, many fits didn’t converge

The plots look OK, but the values don’t make any sense.

x_to_lung <- Xs(myodemodel_to_lung) # make prediction function
loadDLL(x_to_lung)

# get the only the parameters needed for x
pars_to_lung_2 <- rep(0,4) %>% set_names(paste0("k_to_lung_", c("sul", "", "cys", "glu")))

pars <- c(all_pars, pars_to_lung_2)[getParameters(x)]

free_parameters_to_lung_2 <- c("APAPGLU_HLM_CL", "APAPGLU_Km", "APAPSUL_HLM_CL", "APAPSUL_Km", "APAPCYS_HLM_CL", "APAPCYS_Km", 
                             "Ka_apap", 
                             paste0("k_to_lung_", c("sul", "", "cys", "glu"))
                             # "CLrenal_apap", "CLrenal_metc13",  "CLrenal_apap_sul", "CLrenal_apap_cys", "CLrenal_apap_glu", "CLrenal_co2c13"
                             )


mydatalist <- data %>% 
  filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>% 
  filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>% 
  fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name")) %>% 
  select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")


observables_to_lung_2 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
                 apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
                 apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
                 apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_to_lung_2 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)

# free_parameters_to_lung_2 <- c(free_parameters_to_lung_2, scale_parameters_to_lung_2)

# error_model_to_lung_2 <- c(apap = "srel_apap*apap^2 +s0_apap", 
#                  apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu", 
#                  apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul", 
#                  apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_to_lung_2 <- setdiff(getSymbols(error_model_to_lung_2), names(error_model_to_lung_2)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
  cond <- unlist(conditions[i,])[2:3]
  
  trafo <- as.character(pars) %>% set_names(names(pars))
  trafo[names(cond)] <- cond
  trafo[free_parameters_to_lung_2] <- paste0("exp(log", free_parameters_to_lung_2, ")")

  scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_to_lung_2, x = scale_parameters_to_lung_2, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_to_lung_2)
  scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  # scales[1:length(scales)] <- "1"
  
  # errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_to_lung_2, x = error_parameters_to_lung_2, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_to_lung_2)
  # errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  
  trafo <- c(trafo, scales)#, errors)
    
  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p_to_lung_2 <- NULL
for(i in 1:length(p_list)) { p_to_lung_2 <<- p_to_lung_2 + p_list[[i]]}

g_to_lung_2 <- Y(observables_to_lung_2, x)#, parameters = c(free_parameters_to_lung_2, scale_parameters_to_lung_2))  

# err_to_lung_2 <- Y(error_model_to_lung_2, g_to_lung_2)

obj_to_lung_2 <- normL2(mydatalist, (g_to_lung_2*x_to_lung*p_to_lung_2))#, errmodel = err_to_lung_2)


pouter_to_lung_2 <- rep(0, length(getParameters(obj_to_lung_2))) %>% set_names(getParameters(obj_to_lung_2))
pouter_to_lung_2[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
# obj_to_lung_2(pouter_to_lung_2)
# job_to_lung_2 <- runbg({myfit <- mstrust(objfun = obj_to_lung_2, center = pouter_to_lung_2, sd = 2, studyname = "methacetin", cores = 24, fits = 100, iterlim = 150); myfit}, machine = "ruprecht1", filename = "job_to_lung_2", input = global_env_without(c("mypred", "job", "myfit", "myplot")))

# save(job_to_lung_2, file = "job_to_lung_2.rda")

# job_to_lung_2$check()
# myfit_to_lung_2 <- job_to_lung_2$get()$ruprecht1
# save(myfit_to_lung_2, file = "myfit_to_lung_2.rda")
# job_to_lung_2$purge()
load("myfit_to_lung_2.rda")
loadDLL(x_to_lung)
## The following local files were dynamically loaded: methacetin_to_lung.so, methacetin_to_lung_s.so
myfit_to_lung_2 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_to_lung_2 %>% as.parframe())+scale_y_log10()

mypred_to_lung_2 <- (g_to_lung_2*x_to_lung*p_to_lung_2)(mytimes*4, myfit_to_lung_2 %>% as.parframe() %>% as.parvec(35) %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
# myplot_to_lung_2 <- plotCombined(mypred_to_lung_2, mydatalist, name %>% sapply(. %>% str_detect(c(paste0("^",names(observables)), "Agu", "Ali", "Ave")) %>% any))
myplot_to_lung_2 <- plotCombined(mypred_to_lung_2, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_to_lung_2)
# myplot_to_lung_2
##   index    value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1    61 811.8280     FALSE        150         -6.998626     -17.50336
## 2     7 812.1984     FALSE        150         -7.019726     -16.97033
## 3    79 812.3622     FALSE        150         -7.031897     -16.94359
## 4    46 812.3956     FALSE        150         -7.033562     -16.92471
## 5    36 822.9495     FALSE        150         -7.920227     -17.47499
## 6    54 831.7163     FALSE        150         -9.112826     -16.08981
##   logAPAPSUL_HLM_CL logAPAPSUL_Km logAPAPCYS_HLM_CL logAPAPCYS_Km
## 1         -4.647647     -7.067046         -8.005360     -25.35451
## 2         -4.677744     -7.095040         -8.017149     -24.29624
## 3         -4.699630     -7.149380         -8.015090     -24.25044
## 4         -4.703237     -7.161027         -8.014722     -24.18622
## 5         -4.908147     -8.078681         -7.579660     -23.58753
## 6         -5.018467     -8.363523         -7.356687     -23.10982
##   logk_to_lung_sul logKa_apap logk_to_lung_ logk_to_lung_cys
## 1        -2.288579  -8.479607     -4.093847        -2.952564
## 2        -2.275876  -8.451360     -4.094141        -2.997499
## 3        -2.269219  -8.442444     -4.094653        -3.017023
## 4        -2.268317  -8.441194     -4.094752        -3.020282
## 5        -2.031370  -8.290519     -4.068860        -3.637812
## 6        -1.887041  -8.282598     -4.051058        -3.751590
##   logk_to_lung_glu log_scale_apap_Chan1997_14_0
## 1         6.776023                     2.035524
## 2         6.772015                     1.963864
## 3         7.088785                     1.967955
## 4         6.708936                     1.971821
## 5         6.577030                     2.122688
## 6         5.219719                     2.144451
##   log_scale_apap_glu_Chan1997_14_0 log_scale_apap_sul_Chan1997_14_0
## 1                         6.495936                        -2.067167
## 2                         6.514004                        -2.111232
## 3                         6.844040                        -2.125823
## 4                         6.465985                        -2.127919
## 5                         7.263784                        -2.351949
## 6                         7.126559                        -2.320297
##   log_scale_apap_cys_Chan1997_14_0 log_scale_apap_Chiew2010_56_0
## 1                        -4.984626                    -0.1457100
## 2                        -4.992699                    -0.2872532
## 3                        -5.003292                    -0.3473295
## 4                        -5.005029                    -0.3556862
## 5                        -5.625267                    -1.0307136
## 6                        -5.839268                    -1.2249324
##   log_scale_apap_glu_Chiew2010_56_0 log_scale_apap_sul_Chiew2010_56_0
## 1                          8.738710                         -2.787880
## 2                          8.752954                         -2.781700
## 3                          9.079320                         -2.776042
## 4                          8.700771                         -2.775228
## 5                          9.406054                         -2.621662
## 6                          9.203844                         -2.490344
##   log_scale_apap_Critchley2005_14_0 log_scale_apap_glu_Critchley2005_14_0
## 1                          1.914790                              6.775920
## 2                          1.843805                              6.792860
## 3                          1.848134                              7.121887
## 4                          1.852019                              6.743747
## 5                          2.008014                              7.506618
## 6                          2.029117                              7.351299
##   log_scale_apap_sul_Critchley2005_14_0
## 1                             -2.517989
## 2                             -2.561130
## 3                             -2.575354
## 4                             -2.577414
## 5                             -2.786350
## 6                             -2.748679
##   log_scale_apap_cys_Critchley2005_14_0 log_scale_apap_Rawlins1977_0_1
## 1                             -5.196222                      -2.667753
## 2                             -5.205343                      -2.668026
## 3                             -5.216401                      -2.668411
## 4                             -5.218228                      -2.668487
## 5                             -5.860658                      -2.650692
## 6                             -6.081074                      -2.638401
##   log_scale_apap_Rawlins1977_1_0 log_scale_apap_Rawlins1977_2_0
## 1                       2.251798                       1.204153
## 2                       2.148307                       1.137746
## 3                       2.142583                       1.140532
## 4                       2.145176                       1.143994
## 5                       2.125534                       1.249054
## 6                       2.134132                       1.174536
##   log_scale_apap_Rawlins1977_05_0
## 1                        6.104588
## 2                        5.154957
## 3                        4.939859
## 4                        4.913994
## 5                        3.014949
## 6                        2.848154
##                     logAPAPCYS_HLM_CL:  0.000333669227703504
##                         logAPAPCYS_Km:  9.74264737617257e-12
##                     logAPAPGLU_HLM_CL:  0.00091313568705377
##                         logAPAPGLU_Km:  2.50258813575314e-08
##                     logAPAPSUL_HLM_CL:  0.00958412454265644
##                         logAPAPSUL_Km:  0.000852748222445177
##                            logKa_apap:  0.000207660238495109
##                         logk_to_lung_:  0.016674967733492
##                      logk_to_lung_cys:  0.0522056534562965
##                      logk_to_lung_glu:  876.576074026939
##                      logk_to_lung_sul:  0.101410478546347
##          log_scale_apap_Chan1997_14_0:  7.65626070083065
##         log_scale_apap_Chiew2010_56_0:  0.864408357893709
##     log_scale_apap_Critchley2005_14_0:  6.78551630713505
##      log_scale_apap_cys_Chan1997_14_0:  0.00684233640219463
## log_scale_apap_cys_Critchley2005_14_0:  0.00553744269024556
##      log_scale_apap_glu_Chan1997_14_0:  662.443889437685
##     log_scale_apap_glu_Chiew2010_56_0:  6239.84412853667
## log_scale_apap_glu_Critchley2005_14_0:  876.485712326581
##        log_scale_apap_Rawlins1977_0_1:  0.0694080270944431
##       log_scale_apap_Rawlins1977_05_0:  447.908262977161
##        log_scale_apap_Rawlins1977_1_0:  9.50480772104851
##        log_scale_apap_Rawlins1977_2_0:  3.33393329597214
##      log_scale_apap_sul_Chan1997_14_0:  0.126543722411801
##     log_scale_apap_sul_Chiew2010_56_0:  0.0615515731003974
## log_scale_apap_sul_Critchley2005_14_0:  0.0806215896552064

Estimate dosing - some things look OK

Here I try to estimate the dosing. Dunno if this makes much sense, but it might be worth a try.

load("methacetin.rda")

x <- Xs(myodemodel) # make prediction function
loadDLL(x)

# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]

free_parameters_est_input <- c("APAPGLU_HLM_CL",  # Vmax value
                        "APAPGLU_Km",  # Km value
                        "APAPSUL_HLM_CL",  # Vmax value
                        "APAPGLU_Km",  # Km value
                        "APAPCYS_HLM_CL",  # Vmax value
                        "APAPCYS_Km",  # Km value
                        
                        "Ka_apap" ,
                        "D_apap",
                        
                        "CLrenal_apap", "CLrenal_apap_cys", "CLrenal_apap_glu", "CLrenal_apap_sul"
                        #, #"F_apap_sul"    ,
                        # "Kpre_apap", "Kpki_apap", "Kpli_apap",
                        # "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
                        # "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
                        # "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
                        # "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
                        # "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)

fixed_parameters_est_input <- pars[!(names(pars)%in%c(free_parameters_est_input,names(f)[1]))] %>% names

mydatalist <- data %>% 
  filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>% 
  filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>% 
  fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name")) %>% 
  select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")


observables_est_input <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
                    apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
                    apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
                    apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_est_input <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)

# free_parameters_est_input <- c(free_parameters_est_input, scale_parameters_est_input)

# error_model_est_input <- c(apap = "srel_apap*apap^2 +s0_apap", 
#                  apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu", 
#                  apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul", 
#                  apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_est_input <- setdiff(getSymbols(error_model_est_input), names(error_model_est_input)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
  cond <- unlist(conditions[i,])[2:3]
  
  trafo <- as.character(pars) %>% set_names(names(pars))
  trafo[names(cond)] <- cond
  trafo[free_parameters_est_input] <- paste0("exp(log", free_parameters_est_input, ")")
  
  if(cond[1]==0) {
    trafo["Ave_apap"] <- paste0("exp(logAve_apap_", rownames(conditions)[i], ")")  %>% str_replace_all("\\.", "")
    trafo["D_apap"] <- "0"
  } else {
    trafo["D_apap"] <- paste0("exp(logD_apap_", rownames(conditions)[i], ")")  %>% str_replace_all("\\.", "")
    trafo["Ave_apap"] <- 0
  }
  
  
    scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_to_lung_2, x = scale_parameters_to_lung_2, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_to_lung_2)
  # scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  scales[1:length(scales)] <- "1"
  
  # errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_est_input, x = error_parameters_est_input, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_est_input)
  # errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  
  trafo <- c(trafo, scales)#, errors)
  
  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p_est_input <- NULL
for(i in 1:length(p_list)) { p_est_input <<- p_est_input + p_list[[i]]}

g_est_input <- Y(observables_est_input, x)#, parameters = c(free_parameters_est_input, scale_parameters_est_input))  

# err_est_input <- Y(error_model_est_input, g_est_input)

obj_est_input <- normL2(mydatalist, (g_est_input*x*p_est_input))#, errmodel = err_est_input)


pouter_est_input <- rep(0, length(getParameters(obj_est_input))) %>% set_names(getParameters(obj_est_input))
pouter_est_input[names(myfit_CLrenal %>% as.parframe()  %>% as.parvec())] <- myfit_CLrenal %>% as.parframe()  %>% as.parvec()
pouter_est_input <- pouter_est_input[order(names(pouter_est_input))]
pouter_est_input[c(6,11:16)] <- log(c(1,1.4,5.6,1.4,0.5,1,2))

# p_est_input(pouter_est_input)

# job_est_input <- runbg({myfit <- mstrust(objfun = obj_est_input, center = pouter_est_input, sd = 5, studyname = "methacetin", cores = 16, fits = 200); myfit}, machine = "ruprecht2", filename = "job_est_input", input = global_env_without(c("mypred", "job", "myfit", "myplot")))

# save(job_est_input, file = "job_est_input.rda")

# job_est_input$check()

# (g_est_input*x*p_est_input)(times = mytimes, pouter_est_input, deriv = F) %>% plotPrediction(name %in% names(observables))
# load("job_est_input.rda")
# myfit_est_input <- job_est_input$get()$ruprecht2
# save(myfit_est_input, file = "myfit_est_input.rda")
# job_est_input$purge()
load("myfit_est_input.rda")
loadDLL(x)
## The following local files were dynamically loaded: methacetin.so, methacetin_s.so
myfit_est_input %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_est_input %>% as.parframe())+scale_y_log10()

mypred_est_input <- (g_est_input*x*p_est_input)(mytimes*4, myfit_est_input %>% as.parframe() %>% as.parvec(35) %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
# myplot_est_input <- plotCombined(mypred_est_input, mydatalist, name %>% sapply(. %>% str_detect(c(paste0("^",names(observables)), "Agu", "Ali", "Ave")) %>% any))
myplot_est_input <- plotCombined(mypred_est_input, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_est_input)
# myplot_est_input
##   index    value converged iterations logAPAPCYS_HLM_CL logAPAPCYS_Km
## 1    58 1122.299     FALSE        100         -10.48928     -4.957525
## 2   156 1122.300      TRUE         90         -10.49027     -4.961741
## 3    91 1142.864      TRUE         54         -10.48757     -5.127342
## 4   192 1142.864      TRUE         77         -10.48709     -5.124875
## 5   130 1142.864      TRUE         43         -10.48786     -5.128740
## 6   144 1142.864      TRUE         75         -10.48767     -5.127869
##   logAPAPGLU_HLM_CL logAPAPGLU_Km logAPAPSUL_HLM_CL
## 1          7.615842      10.65704         -5.418030
## 2          6.535269       9.57645         -5.418026
## 3         10.713608      13.62289         -5.292063
## 4         18.348451      21.25771         -5.292033
## 5         12.851771      15.76106         -5.292060
## 6          9.264502      12.17378         -5.292044
##   logAve_apap_Rawlins1977_0_1 logCLrenal_apap logCLrenal_apap_cys
## 1                  0.10147172       -6.716203           -5.779922
## 2                  0.10146728       -6.716345           -5.778142
## 3                  0.02988996      -34.131247           -5.648947
## 4                  0.02989775      -32.298981           -5.650080
## 5                  0.02989460      -31.947337           -5.648331
## 6                  0.02989872      -28.053280           -5.648702
##   logCLrenal_apap_glu logCLrenal_apap_sul logD_apap_Chan1997_14_0
## 1           -5.444322           -4.590098               0.8099026
## 2           -5.444325           -4.590114               0.8099012
## 3           -5.331795           -4.452475               0.7546452
## 4           -5.331763           -4.452424               0.7546769
## 5           -5.331793           -4.452482               0.7546558
## 6           -5.331771           -4.452459               0.7546687
##   logD_apap_Chiew2010_56_0 logD_apap_Critchley2005_14_0
## 1                 2.563208                    0.6401644
## 2                 2.563184                    0.6401447
## 3                 2.467657                    0.5631522
## 4                 2.467679                    0.5631632
## 5                 2.467656                    0.5631472
## 6                 2.467666                    0.5631538
##   logD_apap_Rawlins1977_05_0 logD_apap_Rawlins1977_1_0
## 1                 -0.8721721                0.12574694
## 2                 -0.8721339                0.12575422
## 3                 -0.9409273                0.05866417
## 4                 -0.9409427                0.05867302
## 5                 -0.9409016                0.05867707
## 6                 -0.9408988                0.05868636
##   logD_apap_Rawlins1977_2_0 logKa_apap
## 1                 0.8141646  -7.855334
## 2                 0.8141561  -7.855325
## 3                 0.7695203  -7.811703
## 4                 0.7695401  -7.811767
## 5                 0.7695220  -7.811715
## 6                 0.7695283  -7.811742
##            logAPAPCYS_HLM_CL:  2.78331452999271e-05
##                logAPAPCYS_Km:  0.00703030850292055
##            logAPAPGLU_HLM_CL:  2030.10329039433
##                logAPAPGLU_Km:  42490.6189263482
##            logAPAPSUL_HLM_CL:  0.00443587568913188
##  logAve_apap_Rawlins1977_0_1:  1.1067986135081
##              logCLrenal_apap:  0.00121112777686759
##          logCLrenal_apap_cys:  0.0030889577568275
##          logCLrenal_apap_glu:  0.00432076935688847
##          logCLrenal_apap_sul:  0.0101518605119084
##      logD_apap_Chan1997_14_0:  2.24768909406063
##     logD_apap_Chiew2010_56_0:  12.9773760893221
## logD_apap_Critchley2005_14_0:  1.896792740859
##   logD_apap_Rawlins1977_05_0:  0.418042515914904
##    logD_apap_Rawlins1977_1_0:  1.13399516723761
##    logD_apap_Rawlins1977_2_0:  2.25728919295756
##                   logKa_apap:  0.000387678645051781

Do the estimation of _to_lung again with the original data instead of the data with fitted error model - Problem in many fits.

# Multiply rates from "ve" to "lu" with additional free parameters
reactions_to_lung_3 <- reactions
reactions_to_lung_3[reactions_to_lung_3$from %>% str_detect("(1*Ave_apap)"),3] %<>% paste0(" * (k_to_lung_", c("sul", "", "cys", "glu"),")") 

reactions_to_lung_3
el_to_lung_3 <- eqnlist()

for(i in 1:nrow(reactions_to_lung_3)) el_to_lung_3 <- addReaction(el_to_lung_3, reactions_to_lung_3$from[i], reactions_to_lung_3$to[i], reactions_to_lung_3$rate[i], reactions_to_lung_3$description[i])


# Convert to "eqnvec", which is basically a named vector of the ODEs and the names denote the states
f_to_lung_3 <- el_to_lung_3 %>% as.eqnvec()

pars_to_lung_3 <- rep(1,4) %>% set_names(paste0("k_to_lung_", c("sul", "", "cys", "glu")))
parameters <- c(all_pars, pars_to_lung_3)[getParameters(x)] %>% names()
free_parameters_to_lung_3 <- parameters[parameters %>% sapply(. %>% str_detect(c("_CL", "Km", "Ka_apap$", "k_to")) %>% any) ]

# myodemodel_to_lung_3 <- odemodel(f_to_lung_3, modelname = "methacetin_to_lung_3", fixed = setdiff(parameters, c(free_parameters_to_lung_3, "Ave_apap"))) # Include "Ave_apap" because of a bug: You need the sensitivities to at least one initial value
# save(myodemodel_to_lung_3, file = "methacetin_to_lung_3.rda")

load("methacetin_to_lung_3.rda")

x_to_lung_3 <- Xs(myodemodel_to_lung_3) # make prediction function
loadDLL(x_to_lung_3)

# get the only the parameters needed for x

pars <- c(all_pars, pars_to_lung_3)[getParameters(x)]

free_parameters_to_lung_3 <- c("APAPGLU_HLM_CL", "APAPGLU_Km", "APAPSUL_HLM_CL", "APAPSUL_Km", "APAPCYS_HLM_CL", "APAPCYS_Km", 
                             "Ka_apap",
                             paste0("k_to_lung_", c("sul", "", "cys", "glu"))
                             # "CLrenal_apap", "CLrenal_metc13",  "CLrenal_apap_sul", "CLrenal_apap_cys", "CLrenal_apap_glu", "CLrenal_co2c13"
)


mydatalist <- data %>% filter(!is.na(sigma)) %>% select(-n) %>%  as.datalist()
conditions <- mydatalist %>% attr("condition.grid")


observables_to_lung_3 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
                         apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
                         apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
                         apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_to_lung_3 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)

# free_parameters_to_lung_3 <- c(free_parameters_to_lung_3, scale_parameters_to_lung_3)

# error_model_to_lung_3 <- c(apap = "srel_apap*apap^2 +s0_apap", 
#                  apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu", 
#                  apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul", 
#                  apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_to_lung_3 <- setdiff(getSymbols(error_model_to_lung_3), names(error_model_to_lung_3)) %>% set_names(.,.)
p_list <- lapply(1:nrow(conditions), function(i) {
  cond <- unlist(conditions[i,])[2:3]
  
  trafo <- as.character(pars) %>% set_names(names(pars))
  trafo[names(cond)] <- cond
  trafo[free_parameters_to_lung_3] <- paste0("exp(log", free_parameters_to_lung_3, ")")
  
  scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_to_lung_3, x = scale_parameters_to_lung_3, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_to_lung_3)
  # scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  scales[1:length(scales)] <- "1"
  
  # errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_to_lung_3, x = error_parameters_to_lung_3, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_to_lung_3)
  # errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
  
  trafo <- c(trafo, scales)#, errors)
  
  p <- P(trafo, condition=rownames(conditions[i,]))
  return(p)
})
p_to_lung_3 <- NULL
for(i in 1:length(p_list)) { p_to_lung_3 <<- p_to_lung_3 + p_list[[i]]}

g_to_lung_3 <- Y(observables_to_lung_3, x_to_lung_3)#, parameters = c(free_parameters_to_lung_3, scale_parameters_to_lung_3))  

# err_to_lung_3 <- Y(error_model_to_lung_3, g_to_lung_3)

obj_to_lung_3 <- normL2(mydatalist, (g_to_lung_3*x_to_lung_3*p_to_lung_3))#, errmodel = err_to_lung_3)

getParameters(obj_to_lung_3)
 myfit_to_lung %>% as.parframe() %>% as.parvec()

pouter_to_lung_3 <- myfit_to_lung %>% as.parframe() %>% as.parvec()


# job_to_lung_3 <- runbg({myfit <- mstrust(objfun = obj_to_lung_3, center = pouter_to_lung_3, sd = 4, studyname = "methacetin", cores = 24, fits = 1000, iterlim = 150); myfit}, machine = "ruprecht1", filename = "job_to_lung_3", input = global_env_without(c("mypred", "job", "myfit", "myplot")))

# save(job_to_lung_3, file = "job_to_lung_3.rda")

# job_to_lung_3$check()
# load("job_to_lung_3.rda")
# myfit_to_lung_3 <- job_to_lung_3$get()$ruprecht1
# save(myfit_to_lung_3, file = "myfit_to_lung_3.rda")
# job_to_lung_3$purge()
load("myfit_to_lung_3.rda")
loadDLL(x_to_lung_3)
## The following local files were dynamically loaded: methacetin_to_lung_3.so, methacetin_to_lung_3_s.so
myfit_to_lung_3 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_to_lung_3 %>% as.parframe())+scale_y_log10()

mypred_to_lung_3 <- (g_to_lung_3*x_to_lung_3*p_to_lung_3)(mytimes*4, myfit_to_lung_3 %>% as.parframe() %>% as.parvec(35) %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
# myplot_to_lung_3 <- plotCombined(mypred_to_lung_3, mydatalist, name %>% sapply(. %>% str_detect(c(paste0("^",names(observables)), "Agu", "Ali", "Ave")) %>% any))
myplot_to_lung_3 <- plotCombined(mypred_to_lung_3, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_to_lung_3)
# myplot_to_lung_3
##   index    value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1   598 1502.666     FALSE        150          5.222239      9.111693
## 2   585 1502.666     FALSE        150          6.010086      9.899766
## 3   230 1502.691     FALSE        150          7.145232     11.034067
## 4    74 1502.694     FALSE        150          3.539114      7.429065
## 5   901 1502.698     FALSE        150          3.013285      6.905417
## 6   852 1502.701     FALSE        150          2.820840      6.710850
##   logAPAPSUL_HLM_CL logAPAPSUL_Km logAPAPCYS_HLM_CL logAPAPCYS_Km
## 1          4.908114      9.176861         -7.698137     -6.722977
## 2          3.700609      7.968713         -7.698307     -6.723640
## 3          1.294309      5.563838         -7.698262     -6.723571
## 4          2.076930      6.344381         -7.698409     -6.723928
## 5          3.797133      8.063931         -7.696645     -6.715827
## 6          5.819994     10.087409         -7.698370     -6.723744
##   logk_to_lung_sul logKa_apap logk_to_lung_ logk_to_lung_cys
## 1       -0.4319789  -7.093422   -0.05287355         2.874983
## 2       -0.4313477  -7.093479   -0.05289277         2.874915
## 3       -0.4328742  -7.093376   -0.05285579         2.874963
## 4       -0.4307210  -7.093524   -0.05290436         2.874860
## 5       -0.4297375  -7.093220   -0.05277499         2.875319
## 6       -0.4306401  -7.093529   -0.05290571         2.874868
##   logk_to_lung_glu
## 1       -0.7626255
## 2       -0.7628978
## 3       -0.7619719
## 4       -0.7632058
## 5       -0.7651345
## 6       -0.7632765
## logAPAPCYS_HLM_CL:  0.000453671805412169
##     logAPAPCYS_Km:  0.00120295170864251
## logAPAPGLU_HLM_CL:  185.348656922024
##     logAPAPGLU_Km:  9060.6175670097
## logAPAPSUL_HLM_CL:  135.383862495601
##     logAPAPSUL_Km:  9670.74882908226
##        logKa_apap:  0.000830550499923049
##     logk_to_lung_:  0.948499944225183
##  logk_to_lung_cys:  17.7251183888039
##  logk_to_lung_glu:  0.466440188442027
##  logk_to_lung_sul:  0.649223078801294
# save(list= ls(), file = "workspace.rda")

Next possibilities

  1. Take out some data, eg. apap to fit the other parameters such that the fit isn’t dominated by apap (not so good)
  2. Estimate the Km and Vmaxes each individually, then use those values as a prior? (maybe that’s not so good, because: what are the values of the other Kms and Vmaxes?)